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ABSTRACT 


State feedback controllers have been researched extensively in the 
last decade and computational tools have been developed for their design. 
In spite of the better system performance that can be obtained with these 
controllers, their application has not been widespread. This is mainly 
due to insufficient practical design experience and to the scarceness 
of general design guidelines in the published literature. Also, in 
many high-performance applications, especially if state estimates 
rather than states are used for feedback, the system is found to be 
excessively sensitive to parameter variations. 

The purpose of this research is to develop a better engineering 
insight into the design process of state feedback controllers and to 
provide a method for the reduction of their sensitivity to parameter 
variations . 

In this work the design procedure of feedback controllers is des- 
cribed and the considerations for the selection of the design para- 
meters are given. The frequency domain properties of single-input 
single-output systems using state feedback controllers are analyzed, 
and desirable phase and gain margin properties are demonstrated. Spe- 
cial consideration is given to the design of controllers for tracking 
systems, especially those that are designed to track polynomial 
commands. 

As an application example, a controller is designed for a track- 
ing telescope. The telescope has a polynomial tracking requirement and 
possesses some special features such as actuator saturation and multi- 
ple measurements, one of which is sampled. The resulting system has a 
tracking performance that compares favorably with a much more compli- 
cated digital aided tracker. 
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The problem of parameter sensitivity reduction is treated by 
considering the variable parameters as random variables. A perform- 
ance index is defined as a weighted sum of the state and control 
covariances that stem from both the random system disturbances and the 
parameter uncertainties. This performance index is minimized numeri- 
cally by adjusting a set of free parameters. 

A computer program implementing this method was developed and is 
applied to the sensitivity reduction of several initially sensitive 
tracking systems. Sensitivity reduction factors of 2-3 are typically 
obtained with modest increases in output rms and control effort. 
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I. INTRODUCTION 


A. BACKGROUND 

State feedback controllers have been researched intensively in the 
last decade. Design methods for these controllers have been developed, 
the most important among them being quadratic synthesis [BRY-1 ] . The 
application of these methods to all but the simplest systems requires 
computer aid and in recent years, efficient computer programs for their 
implementation have become available [BRY-3]. 

In applying these methods, it is observed in many cases that the 
resulting systems have better performance capabilities than those that 
are designed using classical frequency domain techniques, especially 
for multivariable systems [BU-l]. In spite of this, the application of 
state feedback to practical design problems is not yet widespread, mainly 
due to the fact that this method is relatively new and designers are 
therefore not familiar with its potential. 

In classical frequency domain techniques, a considerable body of 
design experience has accumulated over the years. Based on this exper- 
ience, designers can express the system specifications in terms of the 
controller structure and parameters. No comparable experience exists 
for state feedback controllers and it is therefore difficult, at times, 
to relate the system specifications to design parameters such as state 
and control weights. 

Most of the publications on the subject of state feedback control 
either treat its theoretical aspects or describe the results of specific 
applications. Only few publications such as Bull's [BU-l] give some 
general insight into the design process and provide practical guidelines 
for the selection of the design parameters. Designers experienced in 
the classical techniques are therefore reluctant to use these methods. 

In attempting to apply state feedback to practical designs, the 
problems of sensitivity to parameter variations may be encountered. In 
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most applications of state feedback, almost all the states are estimated 
and not measured. High performance systems using state estimates for 
feedback are typically sensitive to parameter variations. 

The subject of sensitivity has been treated extensively in the 
literature, A collection of papers dealing with this subject has been 
published recently by Cruz [CR-l], Most of the publications, however, 
treat various theoretical aspects of the problem. Design papers, in 
which practical methods for sensitivity reduction are described, are 
relatively scarce. Some of the sensitivity reduction methods that are 
applicable to quadratic synthesis designs are reviewed in Chapter IV 
of this work. None of these methods, however, seems to be directly 
applicable to relatively high order multiinput-output systems in which 
several plant parameters may be variable. The designer is therefore 
apt to compromise on a system that has lower performance but also lower 
sensitivity and will thus not be able to take advantage of the full 
potential of state feedback controllers. 

A better engineering insight into the design process of state 
feedback controllers, and a general method for their sensitivity reduc- 
tion may therefore increase the applicability of these controllers to 
realistic systems. In this thesis, an attempt is made to treat these 
two problems. 


B. THESIS OUTLINE 

In Chapter II the design process of feedback controllers and state 
estimators is reviewed. Various methods for determining the state 
feedback and estimator gains are given and the selection of state and 
control weights for quadratic synthesis is described. 

Some frequency domain properties of single-input single-output 
state feedback controllers are derived. These properties show that these 
controllers, especially if they are designed by quadratic synthesis, 
have advantages over other methods of compensation. 

The application of state feedback to systems with tracking require- 
ments is also described in this chapter. The state augmentations that 
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are required for tracking specific time functions are derived. 

Most of the material in this chapter is taken from sources m the 
literature. Some parts, however, such as the infinite gain margin 
property of full state feedback controllers and the determination of 
the state augmentation required for polynomial tracking, are believed 

to be original. 

In Chapter III, the parameter sensitivity of state estimate feed- 
back controllers is examined and illustrated by means of an example. 

The system in this example is a simplified version of the Stanford 
Relativity Satellite which is a satellite-mounted tracking telescope. 

For this system, it is shown that the sensitivity stems from the use 
of state estimates instead of states for feedback. A stability criterion, 
which is expressed as a frequency margin, is developed for this system. 

In Chapter IV, a method is derived for the sensitivity reduction of 
systems represented in state variable form. The method is applicable 
to multivariable systems with several variable parameters, subject to 
arbitrary inputs. It is an extension of the quadratic synthesis method 
and is based on a method developed by Palsson and Whittaker [PA-l] for 
single-input single-output systems. A computer program implementing this 
method was developed and the results of its application to two systems 
are described. These systems are the simplified and full version of 
the Stanford Relativity Satellite [BU-lj. Considerable sensitivity 
reduction is obtained in both cases, at a modest cost in performance 
and control effort. 

In Chapter V, a state feedback controller for a ground-based track- 
ing telescope is designed. This system was selected as an example for 
the application of the design methods of Chapter II. It has some special 
features such as control through a flexible element, several measurements 
polynomial tracking requirements, and nonlinear elements. Several 
variants of the controller and estimator design are investigated and 
their performances are compared. The sensitivity is reduced using the 
methods derived in Chapter IV, and a nonlinear network is introduced 
for' the improvement of the large signal stability. The tracking per- 
formance is compared to that of an aided tracker designed for a similar 
system. 
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Conclusions and a summary are given in Chapter VI. 

C. CONTRIBUTIONS 

The principal contributions of this thesis are as follows. 

C. 1 State Feedback Control 

(1) Investigation of engineering properties of state feedback 
controllers. Determination of state augmentations required for poly- 
nomial tracking. 

(2) Application of a state estimate feedback controller to a 
high performance tracker. 

C.2 Sensitivity Reduction 

(1) Analysis of the sensitivity of state estimate feedback con- 
trollers and derivation of the frequency margin stability criterion. 

(2) Development of a method for the reduction of sensitivity to 
parameter variations, and of a computer program for the implementation 
of this method. 

(3) Application of the method to estimator design, including the re- 
duced sensitivity design of the attitude control of the Stanford Relativity 
Satellite. 
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II. FEEDBACK CONTROLLER DESIGN 


A. INTRODUCTION 

In this chapter, the design of state feedback controllers for 
linear systems is described and some guidelines are given for their 
application to realistic systems. It is an attempt to fill a gap in 
the literature between the theoretical treatments which typically deal 
with the mathematical steps of designing a controller and the application 
papers which typically describe specific results without discussing the 
design steps and the considerations involved in the design procedure. 

Familiarity with the state variable representation of systems 
and its various conomical forms is assumed. This subject is treated 
extensively in various texts [CA-1, SC-l]. 

Only linear time invariant systems are considered. This is not 
as restrictive as it may seem since the design of feedback controllers 
for nonlinear systems is, in general, very complicated and therefore 
in many cases an open loop controller is designed for the nominal 
trajectory and a feedback controller for the perturbations about this 
trajectory [BRY-2]. The perturbation equations are linear. They may 
not have constant coefficients but in many cases it is preferred to use 
an average value of the coefficients and design a constant gain con- 
troller instead of the much more complex time— varying controller. If 
the range of variation of the coefficients is large and the variation suf 
ficiently slow, a switchable gain controller may be used. 

In the following section, B> regulator design procedures are des- 
cribed, using both pole placement and quadratic synthesis. In the 
quadratic synthesis procedure, the judgement of the designer is involved 
mainly in the selection of state and control weights. Some methods 
are given for this selection based on system requirements. 
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In Section C, some properties of systems using state feedback 
are examined. There seems to be a notion among control designers that 
state feedback may have an advantage over classical design methods for 
multi-input multi-output systems but that it offers no such advantage 
for single-input single-output (SISO) systems. The properties described 
in this Section show that in many cases, state feedback may have con- 
siderable advantages over classical design techniques even for SISO 
systems. 

In Section D, the more complex problems of tracking and nonzero 
mean disturbances are treated and multivariable integral control is 
introduced. 


B. REGULATOR DESIGN 


B. 1 Description of the Design Procedure 

The basic controller for a linear time invariant system is the 
regulator. It is designed to keep the states of the system in an accept- 
able vicinity of zero with an acceptable amount of control, while the 
system is subject to random zero-mean disturbances and the measurements 
are contaminated by random zero-mean noise. 

The system is described by 


x = Fx + Gu + Tw 
y = Hx + v 

where 

x = state vector (n X 1) 

F = system matrix (n X n) 
u = control vector (m X 1) 

G = control distribution matrix (n X m) 
w = state disturbance (q X 1) 

T = disturbance distribution matrix (n X q) 
y = output vector (pxi ) 
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H = output matrix (p X n) 
v = measurement noise (p X 1). 

The regulator design procedure consists of two steps: (1) design 

of a state feedback controller, and (2) design of a state estimator. 

These two design steps are performed independently. 

The state feedback controller has the form 

u = -Gx . 

In its design, all the states are assumed to be available for feedback. 
The design consists in selecting the feedback gain matrix C(m X n) , 

Two principal methods are available for this selection: (1) pole place- 

ment; (2) quadratic synthesis. 

In the first method, the system requirements, such as rise time, 
overshoot, damping ratio, phase margin, etc., are translated into desir- 
able locations of the closed loop eigenvalues and the feedback gains re- 
quired to obtain these eigenvalues are found. Only for SISO systems are 
the gains determined uniquely by the eigenvalues. For multivariable 
systems, additional constraints have to be imposed. Methods for select- 
ing the gains for given eigenvalues are described in Section B.2. 

In the second method, the gains are determined so as to minimize 

a performance index of the form 

00 

f T T 

J = \ (x Ax + u Bu)dt 

*^o 

where the weighting matrices A and B are determined by the designer. 
The methods for selection of the weighting matrices and for determining 
the gains for given weighting matrices are described in Section B.3. 

The estimator generates estimates of the hon-measured states which 
are then used for feedback in the controller instead of the real states. 
In many cases, especially when the system is subject to noise, all the 


- 7 - 



states required for feedback, and not just the non-measured ones, are 
obtained from the estimator. The structure of the estimator and the 
methods for its design are described in Section B.4. 

The engineering properties of systems designed by this method are 
described in Section C. The regulator as described above will keep 
the states of the system in the vicinity of zero when the system is 
subject to random, zero mean disturbances. If the disturbances have a 
constant or slowly time-varying component, or if the output of the 
system is required to follow a command input, regulator type controllers 
may not be able to prevent output errors. State augmentation techniques 
as described in Section D may in some cases extend the regulator design 
procedures to such cases. 

B. 2 Feedback Gain Determination by Pole Placement 

In general, the system specifications will call for either a speci- 
fied time response (rise time, overshoot) to inputs (or recovery from non- 
zero initial conditions) or for a maximum rms level of output and control 
When the system is subject to random disturbances and sensor noise of a 
known intensity. 

Typically, pole placement will be used in the first case since it 
is difficult to correlate the closed loop eigenvalue location with the 
resulting rms value in the outputs and the controls. For single 
input systems, several methods exist for the determination of the gains. 

If the system is in its controllable canonical form [CA-l] the feedback 
gains are simply the differences between the coefficients of the open 
loop characteristic equation and the required closed loop characteristic 
equation [CA-l], If the system matrix has arbitrary form the gain matrix 
C is obtained from the equation 
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c 


( 2 . 2 ) 


= [a - a] T P 

where a - a is the vector of the differences between the coefficients 

of the open loop and the required closed loop characteristic equations t 

P is the transformation matrix from the given form into the controllable 

canonical form [x = Px], Ackerman [AC-l] has provided a formula for 
c 

finding the feedback gains without calculating the open loop character- 
istic polynomial. It is 

C = q T a(F) 
n 

where 

q^ = [0, 0, 1 j e" 1 , 

F = the system dynamic matrix 

a(s) = the desired closed loop characteristic polynomial, 

„ -| 

<*> = [G, FG, . .., F G], is the controllability matrix. 

Some algorithms for gain determination in multi-input system have 
also been developed. Anderson & Luenberger [AN-1 ] transform the system 
into a canonical form that is composed of blocks in the controller 
canonical form on the diagonal and zero blocks above it. The required 
closed loop eigenvalues can be obtained by adjusting the coefficients 
of each block separately by the corresponding control. The arbitrariness 
in the determination of the gains is resolved by engineering considera- 
tions where the designer selects the controls that he prefers to apply. 
Other controls are applied only if the system is uncontrollable by those 
that are preferred. 

Gopinath [GO-l] describes an algorithm for determining the gains 
in which the uniqueness is achieved by restraining the C matrix to 
rank 1. It then has the form 
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c = 


a„ 


a 


[ c -^ j • • • > ] 


where a to a represent the relative contributions of the different 
1 m 

controls and are determined by the designer. For multiple controls^ the 
designer may find it difficult to translate his physical insight into the 
desirable relative contribution of each control. 


Gain determination by pole placement has some basic deficiencies; 
(1) It is convenient to define time response criteria in terms of 
second, or at most, third-order models. The determination of the roots 
of a higher order model that satisfies given criteria is more difficult. 
For a high order system, therefore, there are two possibilities; (a) 

To place all the eigenvalues of the system in the region of the required 
eigenvalues (maybe by assigning double or triple eigenvalues), accepting 
the resulting uncertainty in the time response. (b) To assign all the 
eigenvalues except two or three that are required for the determined 
time response to regions in the s-plane that are far from the imaginary 
axis and thus give faster time responses. This assures that the time 
response will be dominated by the two or three slow eigenvalues. This 
possibility will require higher gains and therefore may saturate the 
controller for even small deflection of the states. (2) In multi-input 
systems, the additional criteria imposed in order to obtain unique gains 
may result in undesirable solutions. For high order systems, pole 
placement therefore does not exploit fully the possibilities of state 
feedback controllers. In general, better designs can be obtained using 
quadratic synthesis. 

B. 3 Feedback Gain Determination by Quadratic Synthesis 

In this method the feedback gains are selected so as to minimize 
a quadratic performance index (PI). For systems which are required to 
recover from non-zero initial conditions and are not subject to further 
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disturbances, the PI is 


J = ^ (x T Ax + u T Bu)dt . (2.3) 

to 

For systems which are subject to random disturbances, the PI is 

t f 

J = lim ■ 1 - - \ (x T Ax + u T Bu)dt . (2.4) 

V" l 

o 

The second case is the more realistic one. A physical interpretation 
of the PI for this case can be obtained by rewriting Eq. (2.4) as 

J m tr [AX(od) ] + tr[BU(oo)] (2.5) 

where 

X(ao) = the steady state state covariance 
U(<=) = the steady state control covariance 
tr = the trace operator. 

The steady state covariance matrices exists since the system is 
assumed to be stable. The PI of (2.5) is thus a weighted sum of the 
state and control covariances where the weighting matrices A and B are 
selected by the designer. 

The quadratic synthesis design procedure consists of two steps: 

(1) selection of the weighting matrices A and B; (2) determination of 
the gains so as to minimize the PI for the selected weighting matrices. 

These steps will now be described. 

( 1 ) Selection of the weighting matrices . The selection of the 

A and B matrices constitutes the actual design in the quad- 
ratic synthesis method since the determination of the gains 
after these matrices are selected is a purely computational 
step. Unfortunately, no set rules exist for the selection of 
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these matrices so as to obtain desired time responses or 
closed loop eigenvalue locations. The selection is, in many 
cases, a matter of designer experience. One of the reasons 
for the scant use of state feedback in actual systems is 
the lack of a sufficient body of this experience among working 
designers. 

In general, the design is iterative, with the weighting matrices 
modified at each step so as to approach the required perform- 
ance. 


Some general guidelines for the initial selection and subsequent 
modification of these matrices can, however, be given. The 
most convenient rule for the initial selection of the weighting 
matrices is probably Bryson's Rule [BRY-l], It is 


a. . 
11 


1 


(X-l ) 

i max 



b. . 
11 


(Ui )‘ 

1 i*ax 


= °> 1 ^ j ; 

-*• J 

h i;i = °< 1 ^ J; 


( 2 . 11 ) 


where xi and u, are the maximum permissible values of 

- L max ‘max 

the respective states and controls. For systems subject to 
random disturbances, they are the rms values of the permissible 
uncertainties in the states and controls. 

This method was used, among others, by Gupta and Bryson [GUP-1] 
and by Bull [BU-l] for the design of controllers. In general, 
it is recommended to weight in the first iteration only the 
outputs and the controls. In most cases, some changes in the 
weights are required after their initial determination. 

For multivariable systems, no general rules can be given for 
the effect of these changes and in general, the selection 
has to be made by trial and error. In some cases, it is 
possible to decouple approximately the system and determine 
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the weights or even design the controllers for the decoupled 
systems. The performance of the overall coupled system has 
to be verified. This method is used by Gupta and Bryson 
[GUP-1]. 

For single input systems, the effect of weight changes on the 
closed loop eigenvalues can be determined by the root square 
locus method. This method is described in detail in Appendix 
c9. It is based on one of the formulations of the quadratic 
synthesis method due to Rynaski and Whitbeck [RY-l], Accord- 
ing to this formulation, a system designed by quadratic syn- 
thesis satisfies the equation 

[B + Y T (-s)AY(s) ]u(s) = 0 , (2.6) 

where Y(s) = (si - F) *G. 

For SISO systems, this equation has 2n roots that are sym- 
metric about the imaginary axis. For a given A and B (a 
scalar for SISO systems), the left half plane roots of this 
equation are the closed loop eigenvalues of the system. 

The effect of weight changes on the closed loop eigenvalue 
locations is determined by fixing the control weight and all 
the state weight except one. A root locus is constructed as 
a function of this weight. If all the weights are varied in 
turn with different values assigned to the remaining weights, 
a grid is obtained that shows the closed loop eigenvalue 
locations as a function of the state weights. The method 
becomes cumbersome if more than two or three states are 
weighted and is totally impractical for multi-input systems 
because of the (m-l)n extraneous roots that are introduced 
where 

m is the number of controls, 

n is the state dimension. 

A design example using this method is presented in Chapter V. 
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The advantage of using the root square locus method instead 
of arbitrary pole placement is that the non— dominant eigen- 
values are not positioned arbitrarily but are located so as 
to optimize the performance according to a PI such as given 
by Eq. (2.3). The dominant eigenvalues have to be determined 
in both cases according to the system specifications. 

Another method that is applicable to SISO systems and that 
does not require the separate examination of the weighting of 
each state is the model performance index method. In this 
method, the state and control weights are determined so that 
the closed loop response of the controlled system approaches 
that of a system having prescribed eigenvalues. It is essen- 
tially an eigenvalue placement procedure in which only the 
dominant eigenvalues and those that are required for the 
cancellation of the zeros are assigned. All others are re- 
moved to regions of non-dominancy in a predictable way. The 
method is based on Eq. (2.6). Several versions of its im- 
plementation are described in the literature [SC-1, RE-l], 

The version presented here is an extension of the method given 
by Anderson & Moore [AN-2, Sec. 5.4]. 

Given the open loop transfer function of the system 

m 

k H (s + z ) 

y(s) _ i=l _ kN(s) (2.7) 

u(s) n D(s) ’ 

yj (s + P i ) 

i=l 

The steps for the determination of A are (b , the scalar con- 
trol weight, is assumed to be unity without loss of generality) 

(a) Define the characteristic equation of the model as 

d m < s > = s r + rj> r 1 + ••• r r • 

The degree of the model must be such that r<n-m-l. 

(b) Define a polynomial 

-14- 



P(s) 


= D (s)N(s) . 
m 

(c) Find a vector d such that 

d T (sI - F) _1 g = (2.8) 

or T 

d ad j (si - F)g = PCs). 

The restriction on the dimension of D (s) is now clear 

T M 

since the dimension of d adj(sl - F)g is less than n - 1. 

(d) Define 

A = dd T . (2.9) 

The results of this procedure can be evaluated using Eq. (2.6) 
(with scalar b). Substituting A from Eq. (2.9) into (2.6), 
it is transformed into 

b + Y T (-s) dd T Y(s) 

T T -1 T -1 

= b + g (-sI-F ) dd (sI-F) g 

_ P(-s)P(s) 

— D -f r_ • 

D(-s)D(s) 

The closed loop eigenvalues for a given value of b will be 
the left half plane eigenvalues of 

bD(-s)D(s) + P(-s)P(s) . 

By constructing a root locus as a function of b, conclusions 
can be drawn about its influence on the closed loop eigenvalue 
locations: (i) For b (high cost of control), the closed 

loop eigenvalues become the stable open loop eigenvalues and the 
mirror images about the imaginary axis of the unstable open loop 
eigenvalues. (ii) For b -> 0 (low cost of control): r eigenvalues 
tend towards the eigenvalues of D^(s); m eigenvalues tend towards 
the left half plane open loop zeros and the mirror images of the 
right half plane open loop zeros; n - m - r eigenvalues form a 
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Butterworth configuration that recedes from the origin as b 
decreases. The distance of the Butterworth eigenvalues from 
the origin tends to [KW-2, Sec. 3.8]: 


U) 


• ■ 0 ?) 


2 \ 2(n-m-r) 


( 2 . 10 ) 


for low b. 

For sufficiently low b the closed loop response of a system 
with left half plane zeros (minimum phase system) will there- 
fore approach the model response. The residues of the eigen- 
values that approach the zeros are small, and the n - m - r 
eigenvalues of the Butterworth configuration are sufficiently 
removed from the origin so as not to influence the time response 
The value of b which will assure the dominancy of the model 
eigenvalues has to be determined for each system individually. 

It is important to note that for systems with right half plane 
zeros (non-minimum phase systems), it is, in general, not possi- 
ble to approach the model response since the m eigenvalues 
that approach the mirror images of the right half plane zeros 
will, in general, have large residues. This is a fundamental 
deficiency of non-minimum phase systems. 

It is also important to bear in mind that the system will have 
the desirable response only to inputs for which the numerator 
of the transfer function is the same as the numerator of the 
open loop transfer functions y(s)/u(s). This includes 
response to set point changes but does not include recovery 
from nonzero initial conditions or disturbances* The time 
response to these latter inputs will generally contain modes 
corresponding to all the eigenvalues of the system, and there- 
fore, in many cases, will be much slower than the response to 
set point changes . 
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For high order systems, it may be undesirable to attempt the 
removal of the Butterworth roots from the region of dominancy, 
since this may require excessive amounts of control (recall 
Sec. B.2). In this case, a zero order model may be used 
which is equivalent to weighting the output only. Obviously 
this is the same as using the root square locus method with 
output weighting only. The time response for low values of 
b will, in this case, be dominated by the n - m ordered 
Butterworth configuration, with ra roots approaching the 
m open loop zeros. The time response of Butterworth config— 
urations is described in many control texts (e.g., KW-2, Sec. 
3.8). The relative weighting of the output and the control 
for this case may be determined by Eq. (2.10). 

There is, at present, no apparent generalization of the model 
PI method to multivariable systems. Anderson and Moore [AN-2] 
suggest a generalization to multi-input systems: by deter- 

mining the ratios of the m components of the control vector, 
an equivalent single input control is created and the method 
may be applied. This is the same constraint that is given 
by Gopinath [GO-l]. The resulting system, however, will not 
be optimal in the sense that the same closed loop behavior can 
be obtained with lower values of the control. In some cases, 
it may be advantageous to sacrifice this optimality in order 
to have a simpler design method , 

The guidelines for the selection of the weighting matrices are 
summarized below. 

(a) Weight only the outputs and the controls for the initial 
iteration. Note that the weight of one of the outputs 
or controls can be set arbitrarily since it is only the 
ratio of the state-to-control weights that is important. 
The initial determination of the weights is most con- 
veniently done by Bryson's rule. 
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(b) For multivariable systems, little can be said about the 
effect of weighting changes except that, in general, in- 
creasing the ratio of output-to-control weighting will 
increase the bandwidth and increasing the weighting on 
states that are derivatives of the outputs will increase 
the damping of the roots, 

(c) For single input systems, root square locus can be used 
for determining the effect of changing the weights of 
the output and perhaps one or two more states. If more 
states have to be weighted, the model PI can be used. 

It is clear from this section that since the weight selection 
is an iterative procedure, it is important to have a fast 
computer program that repeatedly determines the gains and 
eigenvalues for many values of the weights. Such programs 
have only become available in the last few years. 

(2) Feedback gain determination , 'the optimal feedback gain matrix 
that minimizes a PI as given in Eq. (2.3) o- (2.4) is [BRY-2] 

C° = b" 1 G T S° (2.12) 

where S° is the positive definite steady state solution of 
the Ricatti matrix equation: 

S = -SF - F T S + SGB _1 G T S - A . (2.13) 

The optimal feedback is always a full state feedback. 

There exist several computational methods for calculating 
the S° matrix. These methods are described by Bryson & Ho 
[BRY-2], and Bull [BU-l], It seems that currently the most 
efficient method for time invariant systems is the eigenvector 
decomposition method, which has been implemented as a computer 
program (OPTSYS) by Bryson and Hall [BRY-3], This program 
determines the optimal feedback gain matrix and the resulting 
closed loop eigenvalues and eigenvectors from given F, G, A 
and B matrices. Since it is a very fast and inexpensive 
program, it is convenient for determining the state weights 
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iteratively. In general, a satisfactory system will be 
obtained after four to five iterations [BU-l]. 

For low order SISO systems, the gains can also be found by 
hand calculation using an addition to Eq. (2.6) [BRY-1 ] , 
viz. , 

[B+Y T (-s)AY(s)] = [I+C(-sI-F) -1 G] T B(I + C(sI-F) -1 G]. 

For SISO systems, the unique value of C is found by com- 
paring the coefficients of the polynomial that is formed from 
the left half plane roots of the left side to those of 
1 + C(sI-F) _1 G. 

B. 4 Design of Linear Estimators 

As explained in Section B.l, the design of a state feedback controller 
is composed of two distinct steps: (a) controller design, and (b) state 
estimator design. In this section, the design methods for the state estim- 
ator will be described. 

A full state estimator, is a model of the system, the output of 
which is compared to the system output. The difference between the 
outputs is fed back to the estimator through a gain matrix K. For 
the system (repeated from Eq. 2.l) 

x = Fx + Gu + Fw 

y = Hx + v , 

the estimator has the form 

• • ■ 

x = Fx + Gu + K( y - Hx) , (2.1 7-) 

where x is the state estimate^ and K is the estimator gain matrix. 

From Eq. (2.5) it can be seen that the governing equation for the estimate 

^ A . 

error x = x - x is 

* 

x = (F - KH)x+ Fw - Kv. (2.18) 
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T he rate of decay of this error is therefore determined by the n eigen- 
values of F - KH which can be determined arbitrarily by the selection 
of the n X p matrix K, whenever [F, H] is observable. p is the 
number of measurements. Only for p = 1 will the eigenvalues determine 
the components of K uniquely. 

The design of the full state estimator consists of the selection of 
the gain matrix K. Various criteria exist for the determination of K. 
The principal ones are: noise filtration, rate of error decay, and 

sensitivity to plant parameter variations. If the main criterion is 
noise filtration, and if both the process and measurement noise are 
white and Gaussian, the gains can be selected by minimizing a perform- 
ance index of the form [BU-l] 

J 

e 

The estimator that is obtained by this method is a Kalman filter. This 
is an optimal estimator since it generates the least squares estimate of 
x. The Kalman filter gain matrix K° is [BRY-2] 

K° = -PVR- 1 (2.20) 

where P° is the steady state positive definite solution of the 
equation 

P = FP + PF T + Q - PH T R“ 1 HP. (2.21) 

P is the covariance matrix of the estimation error Sf. P° can be 
computed by different methods which are similar to those used for solving 
the Ricatti equation (2.13). Here, too, the most efficient method 
seems to be eigenvector decomposition. The program OPTSYS contains an 
option for evaluating Kalman filter gains. The estimator gains that are 
obtained by this method depend on the relative magnitude of the process 
and measurement noises. Two extreme cases will now be considered. 


lim — — 

t „-*» t -t 
f f o 


T -1 T -1 
(w Q w + v R v)dt. 
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0 


o 

(a) Q = 0; no process noise. In this case, for stable plants, P = 
and K° = 0 are obtained from Eq. (2.21). Since the system is 
not disturbed, the estimate is exact and the noisy error feedback 
can only contaminate it. If the initial state of the plant is 
not identical to that of the estimator, the error decay rate is 
determined by the plant dynamics solely. 

(b) R = 0; no measurement noise. In this case, K°-> » ; i.e,, all 

the eigenvalues of F - K°H tend to - « , which means that the 
estimates are obtained by differentiating the output n - 1 
times. This differentiation can be performed since the output 
contains no noise whatsoever. 

From these extreme cases, it can be seen that as the measurement 
noises decrease relative to the process noises, the estimate error 
eigenvalues of the Kalman filter become faster. If the dominant estimate 
error eigenvalues are much faster than the dominant controller eigen- 
values, noise criteria cease to be important for the determination of the 
estimator gains. In this case the filtering action provided by the 
estimator is not better than that provided by the controller itself. The 
estimator gains can, in this case, be determined by other criteria such 
as error decay rate or parameter sensitivity. 

If error decay rate is the principal criterion, the estimator gains 
may be determined by pole placement. Methods similar to those described 
in Section B.l may be used for this placement. However, if a computa- 
tional method for determining the Kalman filter gains is available, it 
may be more convenient to determine the gains that are required for the 
placement by using this method with artificial noises. The noises in 
this case are merely a computational tool for the determination of the 
gains... 

For single output systems, the influence of the noises on the eigen- 
value locations can be determined by a root square locus procedure. In 
this case the closed loop eigenvalues are the left half plane roots of 

T 

R + z(s)QZ (-s) = 0 ; 
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where 

Z (s) = H(sl - F) _1 r 

R = the measurement noise covariance 
Q .== the process noise covariance. 

For multioutput systems, the eigenvalue locations as a function of 
the noises have to be determinated by trial and error. The method of 
artificial noises can also be used when the eigenvalues obtained from 
the Kalman filter gains are inconveniently located. An example of such 
a problem is given by Bull [BU-lj. In this case the Kalman filter had 
very low damping and could easily become unstable if the filter gains 
were not implemented precisely. The roots obtained by using artificial 
noises are more conveniently located. The resulting system is not 
optimal for the real noises but its rms values are only slightly higher 
than those of the optimal systems. The property of low sensitivity 
of the noise changes to estimator gains has been observed in several 
systems but its generality has not been established. The determination 
of estimator gains from sensitivity considerations is described in 
Chapter IV. 

For a system with p measurements, a full state estimator is 
actually not required for the estimation of all the states. Luenberger 
[LU-2 ] has developed the concept of the reduced order observer which has 
n - p states and which, together with the p measurements, provides an 
estimate of all the n plant states. A method for designing the observer 
by transformation into a canonical form is described in the same paper. 
Gopinath [GO-1 ] developed a design method which does not require such a 
transformation. This method is described in Appendix E. An example of 
the use of such an observer is given in Chapter V. 

If the measurements are noisy, the use of a reduced order observer 
instead of a full state estimator results in larger output and control 
noises. This is so because the noisy measurements are fed into the 
controller directly without being filtered by the estimator. 
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B.5 Performance Index Evaluation 

For a system of the form of Eq. (2.1), which is subject to the 
white Gaussian disturbance vector, w, and for which the measurement of 
the output is contaminated by the white Gaussian noise vector, v, 
the optimal controller is obtained by 

(a) using the optimal feedback gain C° from Eq. (2.12), with 
the feedback obtained from an estimator. Note that the 
optimal feedback is independent of the process and measure- 
ment noise. 

(b) using the optimal estimator gain K° from Eq. (2.20). Note 
that K° is independent of the weighting matrices A and B. 

This is the certainty equivalence principle [BRY-2]. The system is 
optimal in the sense that it minimizes the PI of Eq. (2.5) 

J = tr (AX) + tr(BU). 

It is of interest in many cases to evaluate the consequences of using 
nonoptimal K or C matrices. This can be done by comparing the result- 
ing PI. 

The expressions for X and U that are required for evaluating J 
are given below. If K = K , we have 

X = X + P , (2.23) 

where 

X 

A 

X 
P 

Equation (2.23) is valid because for an optimal estimator 

E[x(t)x T (t)] = 0 . 


lim E[x(t)x T (t)] is the steady state covariance matrix of 
“ _>00 the state 

lim E [x(t )x T (t ) ] is the steady state covariance matrix of 
the estimate 


lim E[x(tjx T (t)] is the steady state covariance matrix of the 
t"* 50 estimation error. 
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The governing equations for x and P are [BRY-2] 

(F - GC)& + X(F - GC) T + KRK T = 0 (2.24) 

and m m 

(F - KH)P + P(F - KH) + KRK + Q = 0. (2.25) 

Note that Eqs. (2.24) and (2.25) are valid for both optimal and non- 
optimal systems but for systems with nonoptimal filters, Eq. (2.23) is 
not valid. For systems with an optimal filter^ the covariance of the 
estimate error may also be determined by Eq. (2.2l). 

For systems with nonoptimal filters, the state covariance is given 
by 

X = X + P + lim E[3((t)x T (t) ] . (2.26) 

In this case, X has to be calculated by using the augmented system 




where 




rp — 

X E(xx ) 

X [2n X 2n] 

A 

_E(xx ) P - 


(2.27) 


(2.28) 
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[2n X q ] 




[2n X p] . 


X is determined from the upper left n X n partition of X A> The 
equation for U is 


A T 

u = cxc . 


(2.29) 


This equation is always valid but if Eq. (2.28) is used to find X, 

X may not be available. In this case 

U = CXC? , (2.30) 

AAA’ 


where 

C = [C, -C ] m X 2n. 

A 

C. ENGINEERING PROPERTIES OF STATE FEEDBACK 
“ CONTROLLERS 

In this section, some frequency domain properties of systems with 
state and state estimate feedback controllers are examined. These 
properties make it possible to compare the results of this design method 
to those of classical frequency domain techniques. This comparison is 
most meaningful for SISO systems, which are the principal domain of 
application of the classical techniques. The eigenvalue separation 
property of state estimate feedback controllers is used in this section. 
This property is derived in most modern control texts [AN-2], It is 
apparent from the state equation for the augmented system of Eq. (2,27), 
which is 
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V 


(2.31) 


x 


X 


F-GC GC 


0 F-KH 


+ 


G 


0 


u + 
o 


U r , 


w + 


0 


K 


where u is an external input that is applied to both the plant and 
o 

the estimator. 


From Eq. (2.31), the characteristic equation of the augmented system 
is 

A(s) = det (si - F + KH) X det(sl - F + GC). (2.32) 

This characteristic equation contains two separate groups of eigenvalues. 

(a) n eigenvalues of the controller (F - GC) that are the same as those 
obtained with feedback from the states instead of from the estimates. 

(b) n eigenvalues of the estimate error system (F - KH). 

It is because of this separation of the eigenvalues that the con- 
troller can be designed as if full state feedback were available and 
the estimator can then be designed separately. Note that this property 
does not depend on the method by which the estimator and feedback gains 
are obtained, but only on the estimator having the structure of Eq. (2.17) 


C.l Transfer Matrices and Tracking Properties 

The transfer matrices from the various inputs to the outputs y 
are found from Eq. (2.31). Defining 


N = si - F + GC 
M = si - F + KH , 

the system output is obtained as 

y = HN _l Gu + HN _1 (sI - F + GC + KH)M _1 rw 
+ HN -1 GCIvf L Kv . 


(2.33) 


From this equation it can be seen that for inputs applied to both the 


- 26 - 


plant and the estimator, the characteristic equation of the transfer 
matrix contains the eigenvalues of the closed loop system only, 

(F - GC), i.e., for these inputs, the system behaves as if the feed- 
back gains were obtained from the states and not from their estimates. 

This is a remarkable property of state estimate feedback design. 

It can also be seen from Eq. (2.36) that the estimate error 5? is not 
excited by the input u q . For all other inputs (disturbances and 
measurement noises), the characteristic equations of the transfer 
matrices contain both the controller eigenvalues and the estimate 
error eigenvalues. 

These different transfer matrices are important for tracker design, 
which is discussed in Section D. 

C.2 Equivalent Compensators 

For SISO systems, it is sometimes of interest to compare the, struc- 
ture of the state estimate feedback controller to that of compensating 
networks obtained by classical design techniques. For this purpose, the 
equivalent compensating networks of the state estimate feedback controller 
are determined. Block diagram representations of the state estimate 
feedback (SEF) controllers are given in Fig. IT-1. 

Figure Il-la is obtained from Eq. (2.1) and (2.17). Figure Il-lb 
is obtained from Fig. Il-la by defining 


G(s) 

= H(sl - 

F) _1 G 

(2.34) 

H (s) 
u 

= C (si 

- F + KH)" 1 G 

(2.35) 

H (s) 

= C(sl 

- F + KH) _1 K . 

(2.36) 


Figure II-lc is obtained from Fig. Il-lb by defining 

H (s) = [1 + H (s)]" 1 . (2.37) 

c u 

Using two linear algebraic relations [GA-l], 
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The eigenvalues of F - KH do not appear in this transfer function. The 

compensating networks H^Cs) an(i ^(s) and the open loop transfer func_ 

tion G (s) can now be used to compare the design with classical frequency 
o 

domain designs. 

C. 3 Gain and Phase Margin 

Single-input single-output (SISO) systems with full state feedback 
controllers have some desirable gain and phase margin properties. Con- 
sider the SISO system 
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X 


- Fx + c Gu 
o 

u = -Cx 

where C is a full state feedback. c q is a variable scalar (it may be 

the gain of the control amplifier) with nominal value of unity. For 

such a controller 

(a) The system has infinite gain margin for changes in c 

-1 ° 

if the zeros of C(sl - F) G are in the left half plane. 
This result remains valid if state estimate instead of state 

feedback is used. The proof of this result is given in 

Appendix A. The system, however, may still be conditionally 
stable, ie. , a decrease in c^ may cause it to become un- 
stable. 

(b) If the feedback gains are found by quadratic synthesis, the 
inequality 

1 1 - C(jwl - F) -1 g| S 1 (2.42) 

is satisfied [AN-2, Sec. 5.3], In fact, it can be shown 
that whenever (2.42) is satisfied, some A and B matrices 
can be found such that the given gains minimize a quadratic 
performance index of the form (2.3). 

If (2.42) is satisfied, the distance of the polar plot of 

C(jwl - F) 1 G from the -1 + jO point is at least unity. 

Since C(jwl-F) 1 G is the open loop transfer function, the 

following results are obtained [AN-2]: (1) the gain margin 

to changes in C q is infinite; (b) the phase margin is at 

least 60° (This is important for determining the tolerance 

of such systems to time delays); (c) conditionally stable 

systems will remain stable for c > | c 

o o nom ♦ 

As an example, the polar plot of a l/s plant which has an 
optimal controller is shown in Fig. II-2. This system has infinite gain 
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margin but is conditionally stable. 



FIG. II-2 POLAR PLOT OF A l/s 3 PLANT WITH OPTIMAL 
CONTROLLER 


D. TRACKERS AND SYSTEMS SUBJECT TO TIME CORRELATED 

DISTURBANCES 

D.l General 

If a system is subject to time correlated disturbances or if it 
is required to track a given reference output, some modification to the 
regulator structure may be required. In this section, reference outputs 
and disturbances are considered which may be represented as outputs of 
linear systems driven by white noise. They are given by 
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(2.43) 


x 

r 


f x + r w 

r r r r 


y 


r 


H x 
r r 


and 


X D - F D X D + 1 D W D 


w = H d x d . 


(2.44) 


For such reference outputs and disturbances, the response may be improved 
by augmenting the state vector and designing a controller for the aug- 
mented state. The required augmentations and the meaning of the improve- 
ment in the controller are discussed in Section D.2. The design of a 
regulator (no tracking requirements) for a system subject to this type 
of disturbance is described by Johnson and Shelton [JO-2]. 

If the intensities of and w^ are zero, a deterministic track- 

ing problem or a regulator with deterministic disturbance result. Only 
the class of F and F Q matrices that cause sustained reference outputs 
or disturbances are of interest in this case; it is treated in Section 


D.3. 

D.2 Stochastic Tracking and Disturbances 


Consider the system of Eq. (2.1), which is required to track a 
reference output of the form of Eq. (2.43) and is subject to disturbances 
of the form of Eq, (2.44). 

An optimal controller can be found for the augmented system 



X 

■ 

'f 0 THp 


X 


0 0 

0 


G 


• 

x 


0 F 0 


X 


r o 

w 

4* 

0 

u 

r 
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-V 


-° ° V 


-V 


-0 r D J 

-v 


-0- 



(2.45) 


y a = 7 " y r = [H ’ “ V 


x 


x 


+ V . 
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So as to minimize the performance index: 

.00 

j = \ [(y - y/vy - y r > + u T Bu]dt. 

J o 

The solution of this problem is given by Kwakernaak and Sivan [KW-2], 

A partitioned Ricatti matrix is obtained and from it the augmented feed 
back gain matrix can be determined as 



C is the feedback matrix from the states of the unaugmented system that 
minimizes the PI 



. 00 

T T 

^ (y A^y + u Bu)dt . 
o 


The optimal feedback is not changed by the tracking requirements or 

the time correlated disturbance. 

C and C are feedforward matrices from the reference and dis- 
r D 

turbance states. These states do not, of course, exist in reality 
and have to be estimated. The block diagram for this system using 
estimators for the states, the disturbances, and the reference states 
is shown in Fig. II-3. 

In Fig. II-3, estimator No. 2 estimates the states x u which are 
defined by 
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y u = [h, o] I 

C A - [C, C D ] . 


The optimal estimator gains K and can be found separately since the 

reference model and the system are not connected. . 

It is important to recognize that a configuration such as shown 
in Fig. II- 3 is not essential for accommodating disturbances given in 
by Eq. (2.43) or tracking reference outputs given by Eq. (2.44), but is 
is the configuration that will give the minimum of the PI of Eq. (2.46). 
Physically that means that the deviations of the output from the reference 
output will tend to be smaller when feedforward as per Eq. (2.47) is 
used. In many cases, the reference output and the outputs of the system 
are not measurable separately and only the error y g = Y “ y r is meas- 
ured. In that case feedforward from the reference output cannot be used. 

A regulator configuration may be used with the output error fed back to 
the estimator. This is shown in Fig. II-4 (with the disturbance omitted). 



FIG. I 1-4 ESTIMATOR WITH OUTPUT ERROR MEASUREMENT 
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The transfer matrix from y to y for this case will be the same 

r 

as the transfer matrix from the measurement noise v to y (Eq. 2.38). 

y = H(sl - F + GC) _1 GC(sI - F + KH) _1 Ky r , (2.48) 

i.e., in this case, the system has no way of distinguishing between the 
reference output and the measurement noise since it has no information 
about the reference model. In general, the output error for this con- 
figuration will be greater than for the optimal configuration of Fig. 
IT-6, 

An exception is the case of complete reducibility of the reference 
output [RA-1 , KR-1 ] . This is defined as the case in which an error 
state vector e = x - x^ can be defined such that 

e = Fe + Gu (2.49) 

y = He. 
e 

For a completely reducible reference output, the error state therefore 
has the same dynamics as the system. The condition for complete reduci- 
bility is that all the eigenvalues of F^ have to be among those of F. 
A state transformation may be required in order to obtain the form of 
Eq. (2.49). If the reference output is completely reducible, feedback 
from the error states (or their estimates): u = - Ce will minimize 

the PI of Eq. (2.46). The configuration of Fig. II-4 is therefore 
optimal for this case. 

D. 3 Deterministic Tracking and Disturbances 

Xf w and w in Eq. (2.43) and (2.44) vanish, a deterministic 
r D 

tracking problem results. This problem can be solved by state augmenta- 
tion in the same way as described in the previous Section, but the re- 
sulting solutions for most forms of F^ and F^ matrices represent con- 
trols that are optimal for specific reference outputs and disturbances 
only and therefore have little general interest. 
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A class of reference outputs and disturbances that is of general 
interest is the class of polynomials in time. They can be considered 
as outputs of models of the form 

F 

r 

The order of the model is determined by the order of the highest polynomial 
output. This output is the first state of the model. The coefficients 
of the polynomial are determined by the initial conditions. This class 
of reference outputs includes steps, ramps, constant accelerations, etc. 
Systems that are required to track this type of reference output or are 
subject to this type of disturbance generally are required to have zero 
output error at equilibrium. 

The interest in this class of inputs comes from the fact that if a 
system is designed for error free tracking of a polynomial reference 
output, it will also track error free all polynomial reference outputs 
of lower order. It is therefore adapted to an entire class of reference 
outputs and not just for one specific function of time. The same is 
true for polynomial disturbances. 

The consideration of polynomial tracking in conjunction with state 
feedback follows the classical design approach in which the type of a 
system is defined as the order of polynomial that it can track with 
finite error. 

f 

The following points now have to be determined. (a) What are 
the requirements in the structure of the system such that it shall be 
capable of error free tracking of a polynomial? (b) What are the addi- 
tions to state feedback required in the controller in order to obtain 
the error free tracking? 


^ The material in this part is an extension of Ch. 3.7, Ref. KW-2. 


0 ! I i-i 

I 

1 
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The equilibrium condition of a system tracking a polynomial ref- 
erence output is defined as the condition in which the error has vanished. 
In this condition 


x (t) = Fx (t) + Gu (t) 

s s s 

y (t) = Hx (t) = mt r , 

s s 


(2.50) 


where x g (t) is ihe equilibrium solution of Eq. (2.1) for the given 
reference output and m is an arbitrary vector [l X p]. 

Whether such a condition can be reached depends on whether Eq. (2.50) 

r 

can be solved for u (t) for given y = mt . To determine the conditions 

s 

for the existence of u g (t), Eq. (2.50) is rewritten in the form 


si - F 

G 


X 

s 


* 0 

H 

0 


u 


m 




s 


_s r+1 J 


A solution for u and x can be obtained if the matrix 
s s 


sI-F G 

H 0 


is square and has full rank. This condition is fulfilled if (a) the 
number of controls equals the number of outputs, and (b) the system is 
controllable and observable. In this case, 

u (s) = [H(sl - F) -1 G] -1 ~— r . (2.52) 

s r+i 

s 

If p> m (more outputs than controls) , no solution exists for u . If 

■ ■ . ■ . s 

m > p (more controls than outputs) , additional outputs may be defined or 
u g can be optimized for some performance index. Holley & Bryson [HO-l] 
describe such an optimization for a steady reference output. 
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The conditions (a) and (b) above guarantee the existence of 
u (t). However, for this control to reach an equilibrium state which is 
a polynomial of time, one more condition is required. The matrix 

U(s) = [H(sl - F) _1 G] _1 

has to be stable. If this is not so, the equilibrium state of zero 
output error is achieved but the control diverges. The additional condi- 
tion can also be posed as the requirement for the open loop transfer 
matrix 

Y(s) = H(sl - F) _1 G 

to have minimum phase which means that all the roots of the numerator 
polynomial 

det[H ad j (si - F)G] 

are required to be in the left half plane. 

It is important to bear in mind that no real system can track 
indefinitely a polynomial reference output that results in polynomial 
u of order greater than zero, since all controllers are subject to 
saturation. Two cases are therefore of practical importance. (1) The 
tracking of function of time that can be assumed to be composed of poly- 
nomial segments of such length that the controller is not saturated. 

(2) the tracking of polynomials that result in constant control. 

If u exists, an error system can be defined as 

s ■■.■■■ 

e(t) = x(t) - x (t) 
s 

u (t) = u(t) - u (t) (2.53) 

e s 

y (t) = y(t) - y (t) . 

e s 
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Substituting Eq. (2.53) into Eq. (2.1) and using (2.50), the error equa- 
tion is obtained as 


e = Fe + Gu 

e 

y e = He . 


(2.54) 


This is a regulator problem for which a state feedback gain can be found 
in the form 

u = - Ce . (2.55) 

e 

The total control that is required for error free tracking is therefore 
the sum of the error feedback and the equilibrium control 


u (t ) = -Ce(t) + u (t) . (2.56) 

s 

According to Eq. (2.54) the error feedback can be obtained from 
an estimator that has the dynamics of the system and the output of 
which is compared to the error output y^. 

Zero equilibrium output error is assured if the equilibrium control 
is obtained by feeding back integrals of the output error. The number 
of integrations that is required can be determined by finding the number 
A for which the solution of the following equation (which is obtained 
from Eq. 2.52) exists: 


lim s^u = lim s^[H(sI - F) 1 G] _1 m — = const. (2.58) 

s->o s s-»o g r+l 

By assigning a nonzero value to the jth component of m only, 
the number of integrations n^ of the jth output error is found as 

n . = A . . (2.59) 

J J 

If u e = const is imposed in (2.57), the order of polynomials that 
can be tracked with constant control is obtained. 
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lxm s [H(sl - F) G] m -j— ■ = const . (2.58) 

s 

The constant control u is obtained by one integration of the output. 

For SISO systems, £ is the number of roots at the origin (free 
integrations) of the open loop transfer functions. Comparing this condi- 
tion with the condition for complete reducibility (Eq. 2.49) it results 
that for polynomial reference output, complete reducibility is equivalent 
to the possibility of error free tracking with constant control. 

In order to determine the gains from the error integrals, a state 
vector z is defined such that 


z s F. z + G. (y - y ) 

1 i J J r 


The order of z is 


n 

z 


The matrix F [n X n ] is defined as 
i z z 




where 


A-x - 
Ai 


°! %-i 

i J 
i 


Lo! o 


(2.60) 
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The matrix G. [n X q] is defined as 
X z 

'g(A v J) ~ 1 

i/v 

. g(i,j) = 0 

The equilibrium control is therefore obtained as 

u = -C z . 
s o 

The sain matrix C [m X n ] may now be determined by quadratic synthesis. 

o z 

The system is augmented by the state z and the control is selected to 
minimize a PI of the form 

J = ( (x T Ax + z T A q Z + u T Bu)dt . (2.61) 

A different method for determining the integral control gains 
was developed by Holley and Bryson [HO-l], In this method the additional 
roots due to the integral control are determined separately after the 
system is designed. 

The block diagram of the complete tracker is shown in Fig. II-5. 



FIG. II-5 TRACKER WITH INTEGRAL CONTROL 
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In Fig. II-5 the integral control is applied to both the plant 
and the estimator. The application to the estimator is not essential, 
but if it is done the eigenvalue separation between the controller and 
the estimator is preserved. The eigenvalues of the controller are now 
those of the augmented system 


F. 


-GC 


G H 
1 


F-GC 


(2.62) 


The design of systems that are subject to polynomial disturbances is 
done along the same lines. Note, however, that it is rare to have 
polynomial disturbances of order higher than zero (constant disturbances). 

For the sake of completeness, the development for the general case 
is given below. 

It is required to determine: (a) the system structure that will 

permit the output error to remain zero when the system is subject to 
polynomial disturbances; (b) the additions to the feedback required 
in order to keep the error at zero. 


The system is defined by 


k = Fx + Gu + rw (2 .63) 

y = Hx , 

wh.Gr© w ~ m *fc « In *th© stGsdy s"t<rt© 

Hx = 0 . (2.64) 

Using (2.64) and the Laplace transform of (2.63), we obtain 

H(sl - F) 1 (Gu + Tw) = 0. (2.65) 

w 

a can now be determined for a given w if the conditions for the 
w 
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to be non-divergent , the 


existence of Eq. (2.52) are satisfied. For U w 
minimum phase condition also has to be satisfied, is then given by 

u = [H(sl - F)' 1 G]' 1 H(sI - F)" 1 Tw . (2.66) 

w 

U can be obtained, as for polynomial reference output, by feeding back 
w 

integrals of the output error. The number of required integrations rt 
is determined by solving 

lim s rt [H(sI - F) _1 G]“ 1 H(sl - F) -1 r = const. (2.67) 

s 

The controller shown in Fig. II-5 can therefore be used for this 
case without change if the number of integrations is determined as the 
larger of those required to track the reference output and to reject 
the disturbances. If the integral control is used for cancelling dis- 
turbances only, and is applied to the estimator , a steady state error 
in the estimator output results [TA-l], The system output will, however, 
be maintained at zero. 
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Ill 


SENSITIVITY OF CONTROLLERS TO PLANT 


PARAMETER VARIATIONS 


A. INTRODUCTION 

The design procedure of state feedback controllers and some of their 
engineering properties were described in Chapter II. The controller and 
estimator gains depend on the structure of the system matrices and on the 
values of their parameters. If the actual values of these parameters 
differ from those that were assumed in the design, some of the system 
properties described in Section II-C may not be preserved. In most en- 
gineering problems, the system matrices are only an approximate des- 
cription of the actual system structure and the values of their parameters 
are, in many cases, ill defined and may also vary in time. The sensitivity 
of the controller to plant parameter perturbations is therefore of great 
importance to the designer. Some aspects of this problem have been 
treated extensively in the literature, using the various definitions of 
sensitivity given in Section III-B. One aspect that has received exten- 
sive treatment is the sensitivity comparison of open loop and closed loop 
controllers [e.g , KR-2], especially if the loop is closed by full state 
feedback or by full state estimate feedback [KW-l], 

This problem, however, is not especially relevant to the designer. 
The decision to use a feedback controller instead of an open loop con- 
troller is, in general, made for various reasons other than sensitivity. 

The problem is then to compare the sensitivities of various types of 
applicable feedback schemes. There is no general solution to this prob- 
lem and each system has to be examined individually. 

In the context of this work, the sensitivity of state estimate 
feedback controllers is of special interest. It is important to disting- 
uish, in this respect, between state feedback and state estimate feedback 
controllers . While the nominal properties of these types of controllers 
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may be close, their sensitivies are generally totally different. In 
most cases, the state estimate feedback (SEF) controllers will be more 
sensitive but this is not a meaningful comparison since the estimates 
are only used because the states are not measured. Therefore the choice 
is not between state feedback and SEF controllers but between SEF con- 
trollers and other applicable feedback schemes. 

In Section III-C the sensitivity of a SEF controller is analyzed 
and in Section IH-D the sensitivities of such a controller and a classi- 
cal compensation network are compared for a specific case. 


B. MEASURES OF SENSITIVITY 

The functionally most significant measure of sensitivity which is 
also most frequently used in time domain design is trajectory sensitivity. 
It is defined as the vector 


cr(t) 


dx(t) 
d (1 


(3.1) 


where [1 is a scalar variable parameter [KR-3], If there is more than 
one variable parameter, a sensitivity vector is defined for each para- 
meter. 


In order to get a scalar, time independent measure of sensitivity, 
the time average of a weighted square of this expression is used in 
stochastic systems and its time integral in deterministic systems. These 
integrals or average values are then defined as sensitivity performance 
indices [KR-3]. 

For deterministic systems 

00 CO 

J SD = J o- T (t)W<j(t)dt = trjw J o cr(t)u T (t)dtj. (3.2) 

For stochastic systems 
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j = lim — \ E[a T (t)Wa(t) ]dt 

ss t -»» t ) 

1 o 

f tf 

= trlw lim ~r~ \ E[a(t)a T (t) ]dt 

l t -*» t f ) 

n 


(3.3) 


where W is an arbitrary weighting matrix. For a stable stochastic 
system, the approximate value of the performance index is 

J gs = tr |wE [ (t(<» ) (°°) ] j 1 , (3.4) 


where cr(«) is the steady state value of the sensitivity vector. The 
governing equations for o'(t) are given in Chapter IV-A. 

Trajectory sensitivity and the sensitivity performance index (PI) 
related to it are good measures of sensitivity but they are difficult 
to evaluate and therefore a simpler measure is often required, especially 
for preliminary design. One such measure is the eigenvalue sensitivity 
defined as [PO-1, Sec. 32] 


where X is the kth eigenvalue, 
k — 

This measure of sensitivity is not as closely related to the 
system time response as the trajectory sensitivity because: (a) only 

the dominant eigenvalues influence the time response; (b) the time 
response perturbation depends not only on perturbations of the eigenvalues 
but also on those of the eigenvectors. In many cases, however, especially 
when different controllers for the same system are compared, a sufficient 
measure of their relative sensitivity can be obtained by comparing the 
parameter perturbations required to induce instability. For this case, 
eigenvalue sensitivites are adequate. 

A numerical approximation of the eigenvalue sensitivities can be 
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obtained by computing the eigenvalues of the nominal and perturbed dynamic 
matrices. The QR Algorithm by Francis [FR-l] is a very efficient method 
for this computation and it has been incorporated into a variety of 
computer programs, among them the OPTSYS program used for optimal con- 
troller design. 


In frequency domain design, transfer function sensitivity is often 
used. It is defined by Bode [BO-l] as 


s = d In T(sjl) (3.6) 

s d in |i 

where 

T(s,|l) = G(s,|i) 

1 + G(s ,|i) H(s ,|i) 

is the closed loop transfer function. This sensitivity measure is not 
used in this work. 

C. S TATE ESTIMATE FEEDBACK CONTROLLERS WITH 
~ “ PERTURBED PARAMETERS 


Consider the system of Eq. (2.1) and (2.17) with the system matrices 
perturbed as follows: 

F = F Q + 5F 

G = G Q + 5G (3.7) 

H = H q + 8H . 


The dynamic equations are now 



r f +&f 

o 

~(G o +5G)c" 

+ 

r" 

w + 

”0“ 

v + 

G +5G _ 
o 

_K(H q +5H) 

F -G C-KH _ 
o o o 


-0_ 


_K 


- G - 
o 


(3.8) 
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It is assumed here that the errors occur only in the matrices F, G, and 
H over which the designer has no control. No errors are assumed in the 
gain matrices K and C. If such errors exist, their influence can be 
analyzed by replacing SGC by B(GC), and K5H by 5(KH) in Eq. (3.8). 

Using the transformation 


r* 


P 

- 


" — 

X 


I 

0 


X 

5c 


I 

-I 


/S. 

X 

_ — 


_ 

- 


M* - 


Eq. (3.8) is changed to 



F +5F-(G +6G)C (G +SG)C 

o o o 


" - 
X 

4* 

r 

w - 

0 

v + 

G +SG 
o 

5F-5GC-K6H F -KH +5GQ 

L O O 


X 


_r_ 


K_ 


SG 

L J 


(3.9) 


The salient fact that can be observed from comparing Eq. (3.9) 
and (2.36) is that the parameter perturbations couple the state into the 
estimate error equation and therefore destroy the separation of the con- 
troller and the estimate error eigenvalues. As demonstrated in III-D, 
this coupling may lead to instability for even relatively small perturb- 
ations if the original eigenvalues have low damping. 

The characteristic equation of the perturbed system may be found 
by using the expression for the determinant of block matrices [GA-1, Sec 
5] 

= det M det(M 4 - = 

= det M det (M^ - m 2 m 4 1m 3 ^ • 

Defining 

M = si - F q + G q C 

N = si - F q + KH q - 5GC .. 


det 


M 1 M 2 


M M 
3 4 
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The characteristic equation of (3.9) is obtained as 


D(s) 


det 


M-&F+SGC -(G +5G)C 
• o 


-5F+SGC+K5H 


N 


- 0 


(3.10) 


which reduces to 


D(s) 


det N X det [M - (sI-F +KH +G C)N -1 (SF-&GC) 
v o o o 

+ (G o +SG)CN"' 1 K&H] = 0 . 


(3.11) 


The influence of the parameter changes on the system eigenvalues may be 
determined from Eqs. (3.9) or (3.11). Equation (3.11) can be put in the 
form of a sum of terms of which one is the unperturbed characteristic 
equation and the others are the coefficients of the parameter perturba- 
tions, viz., 


D(s) = D (s) + D 1 (s)8|JL 1 + + D r (s)5|l r =0. (3.12) 

From this form the influence of each parameter may be evaluated by root 
locus techniques. Computer programs are available for this evaluation. 

The root loci may also be constructed directly from Eq. (3.9) by perturbing 
the parameters and computing the eigenvalues. 

D. A NUMERICAL EXAMPLE OF THE SENSITIVITY OF 
STATE ESTIMATE FEEDBACK (SEF) CONTROLLERS 

D.l General 

In this section a detailed numerical example is given of a state 
estimate feedback <SEF) controller for which the parameter sensitivity 
problem is particularly severe. This same system is used in Chapter IV 
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for an application of the sensitivity reduction method which is developed 
there. This system was selected because it demonstrates the sensitivity 
problem of the class of systems in which there is an elastic element 
between the controller and the controlled element. Also, it is suffi- 
ciently simple so that the underlying reasons for its sensitivity can 
be perceived. 


D.2 Plant Description 


The plant used in this example is derived from the Stanford 
Relativity Satellite, a controller for which was designed by Bull [BU-l], 
This Satellite, which is shown in Fig. III-l, consists of an outer body 
in which a helium filled dewar is mounted elastically. The dewar con- 
tains a telescope which is also connected to it elastically. The 
attitude of the telescope is controlled by means of two controllers: 
an actuator between the telescope and the dewar, and a thrustor mounted 
on the outer body. The actuator provides the high bandwidth precision 
control for the telescope and the thrustor controls the outer body so 
that the relative attitude between the three bodies remains small. The 
attitude of the telescope is measured. 

A dynamic model of the plant is shown in Fig. III-2. The damping 
constants b^, and b^ are very low and can be neglected in the analysis. 
The spring constants k and k are poorly defined and may vary in time. 

y 

The controller gains, especially for the thrustor, may also vary in time, 
but the sensitivity to these variations is low and they are therefore not 
considered here. In this example, a low frequency approximation of this 
plant is used. For low frequency behavior, it can be assumed that the 
spring force in the spring k-y is cancelled by the control so that 

no net force is applied by the telescope on the dewar. 

The model of the approximated system is shown in Fig. III-3. Its 
governing equations are 


1^0 = kcp + w^ 

I (cp + 9 ) = -kcp + u - w + w 

2 * 1 ^ 

The state vector x is defined as 


(3.13) 


- 51 - 


SOLAR ARRAY 
’(4 PANELS) 


ACQUISITION 
SUN TRACKER 


2 ANTENNAS 
LOCATED £N 
OPPOSITE 
PANELS 


CONTROL NOZZLES 
(6 PAIRS) 


ACQUISITION . 
STAR TRACKER 










QUARTZ BLOCK 
ASSEMBLY 


LIQUID 

HELIUM 


SUN SHIELD 


ENTIRE GIRTH RING 
■REMOVABLE FROM 
DEWAR AS A UNIT 


ELECTRONICS 
'(4 BOXES) 
ON THERMAL 
STAND-OFFS 






^•X\V\wVv 




WM, 


FIG. III-l LAYOUT OF THE STANFORD RELATIVITY SATELLITE 
[ from BU-l] . 
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Using k/l = w 2 (the natural frequency of the oscillatory motion), 

and k/ I]L = 0^2° where I3 - + V* the State re P resentation 

of Eq. (3.13) is obtained. 
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(3.14) 


The disturbances that are acting on this sytem are [BU-l]: (a) d^, 
thrustor noise torque acting on the outer body at an equivalent bandwidth 
of 1 rad/sec. (b) 0^, noise torque between outer and middle body at an 
equivalent bandwidth of 1 rad/sec. 

Their intensity matrix has the form 


2 


* n \2 /d^ 2 ‘ (3.15) 

V + '< r 2) 

The intensities of these disturbances are not well known and in the actual 
computation of Bull [BU-l], the off-diagonal terms were neglected and the 
diagonal terms were adjusted so as to obtain desirable pole locations. The 
intensity matrix is thus considered as a pole placement device. 
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1 






In order to obtain the same estimator roots and sensitivity prop- 
erties, this same intensity matrix is used here. 

The numerical values of the system parameters and the covariance 
matrix elements are thus 


aw 

0) 

i 


2 

0 

2 

0 


2 


= 19.5 sec 

oc “2 
= 25 sec 

2 

= 250 kg m 


Q. 


w 


-14 

1.1X10 0 


0 2.3X10 


■14 


[rad^ sec 


The measurement noise intensity is 

-17 2 

r = 5.5 x 10 rad sec. 


D.3 Contoller and Estimator Des ign 

An optimal controller for the system was designed using the weight - 


ing matrices derived 

from 

Bull [BU-1] 





~12X10 4 0 

0 

0 



0 5 

0 

0 

A 

= 

0 0 

6.5X10 5 

0 



0 0 

0 

6.5X10 


The control gains were found using the OPTSYS program [BRY-3]. The 
optimal estimator gains were found using the same program, with the 
disturbance covariance matrix of Eq. (3.15) and the given measurement 
noise variance. 
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The optimal gains are 


C = [109, 270, 284, 384] 

and 

K T = [18.7, 175, 13.9, -30.5]. 

The eigenvalues are: (a) controller (F - GC) 

- 0. 35 + 5. Oj 

- 0.41 + 0.41 j ; 

(b) estimate error (F - KH) 

- 1.01 + 5 , 11 j 

- 8.33+ 8. 35 j . 

D.4 Sensitivity Root Locus 

Since there is only one significantly varying parameter (the 
spring constant k) , the characteristic equation of the perturbed system 
can be written in the form of Eq. (3.12) as 

D(s ) - D Q (s) + D 1 ( s ) & k • 

The eigenvalues of D Q (s) are those of the optimal controller and 
estimator. D (js) is a sixth order system with the eigenvalues: 

-0.41 ± 0. 41j , -2.21 ± 5 . 26 j , -7.47 ± 7.68j. With these eigenvalues, 

a root locus as a function of Sk was constructed. It is shown in 
Fig. III-4. Only the regions of the root locus which are in the vicinity 
of the nominal values are significant since very large parameter 
perturbations are generally not expected. For this system, the region 
of k < 0 has no meaning. 

The range of stability of the system is from the root locus: 
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cr (sec -1 ) 


FIG. I I 1-4 SENSITIVITY ROOT LOCUS FOR STATE ESTIMATE 

FEEDBACK CONTROLLER 
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= + .26 


jw(sec -1 ) 


0 73 < — < 1.26; or, 18.1 < k < 31.5. (3.17) 

k o 

It is important to note that this is not the permissible operational 
range since the performance will generally deteriorate before the system 

becomes unstable. 

While the root locus gives a good representation of the sensi- 
tivity problem, it does not provide a method for its solution. It is 
not obvious how the system eigenvalues have to be modified m order to 
decrease the sensitivity, especially since the root locus numerator 
eigenvalues are also determined by the gain matrices K and C, and 
therefore will change whenever the system eigenvalues are changed. 

For a SISO system, more insight can be gained into the sensitivity 
problem by considering the transfer function and using frequency 
domain stability criteria. In principle, this representation may also 
suggest the modifications required for reducing the sensitivity but 
in practice the implementation of these modifications is difficult, 

D.5 Frequency Domain Analysis 

From Eq. (2.41) the equivalent open loop transfer function for 
a system using an estimator is 

G o (s) = G c (s) G p (s) , 

where 

G (s) = H(sl - F) G 

P 

is the plant transfer function, and 


. . C ad.j[sl - F + KH]G 

G c (s) = det [si - F + KH + _ GC] 

is the equivalent compensator transfer function. 

For the example, these transfer functions were calculated using 
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the computer program XAGSA [WI-l], They are 


G (s ) 
P 



G (s ) 
c 


4.1 X 10 5 (s + 0.37)(s 2 - 1.58s -f 5,2 2 ) sec " 2 # 
(s 2 + 3.4s + 5.85 2 )(s 2 - 1.65 + 11.8 ) 


G Q (jw), the open loop frequency response, is shown in Fig. HI-5, 

v ' Since the plant is undamped, the frequency response has three zero 

crossings; i.e., three points at which the gain equals unity (points 
A, B, C). The frequency domain stability criterion for a system 
with several zero crossings can be found from its polar plot. The sig- 
nificant parameters for determination of stability are the phase 
angles at the three zero crossings. They are shown in Fig. II 1-5 

Polar plots are sketched in Fig. III-6 for four different cases 
of the angles co and r p with <p < 180° in all cases. The cases 

are: 

(a) cp < 180° ; 

D 

(b) cp B > 180°; 

(c) cPg > 180° ; 

(d) < 180° ; 

Case b is shown in Fig. III-5. Case a is obtained if the gain for 
this system is increased. Cases c and d cannot be obtained by gain 
changes. Case c can be obtained by parameter variation (see below). 
The exact shape of the polar plots is not important for the stability 
determination as long as the quadrants of the zero ci’ossings are pre- 
served . 

From Fig. III-6 it can be seen that only case b is stable. 

The condition of stability for this system is therefore 


cp c < 180° 

ip < 180° 

■C 

cp c > 180° 

T c > 180 ° . 
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FIG. III-5 OPEN LOOP FREQUENCY RESPONSE OF STATE ESTIMATE 
FEEDBACK (SEF) CONTROLLER 
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cp > 180 

B 


(3.18) 


cp c < 180° . 

The sensitivity to variations in can be determined from this 

criterion. 

Figure III-7 is an amplitude frequency plot of the region 
of oo = with part of the phase-frequency plot overlaid. 



FIG. III-7 FREQUENCY RESPONSE IN THE REGION OF RESONANCE 


Since the plant root at w = has no damping, its influence on the 
phase frequency plot consists of the addition of a phase lag of 180° 
at this frequency without modifying the plot at other frequencies 
Also, the shape of the amplitude frequency plot in the vicinity of the 
frequency is determined mostly by this root. Changes in the natural 
frequency of the plant will therefore cause the amplitude plot to move 
relative to the phase plot in Fig. III-7 without changing the shape of 
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the phase plot and with little change in the shape of the amplitude 
plot. It can be seen from the figure that changes in u> 0 will cause 
instability as follows: decrease by 1 rad/sec instability as per 

Fig. III-6a; increase by 0.8 rad/sec -+ instability as per Fig. III-6c. 

The region of stability found by this method is somewhat larger 
than that found by the root locus method (Eq. 3.16). In order to 
decrease the sensitivity, the frequency margins 5^ and 5w u have to 
be increased. This can be done by modifying the compensator roots so 
that the phase slopes in the region of the zero crossing B and C are 
decreased. However, decreasing the slope in the vicinity of B will 
also decrease the phase margin at the first zero crossing A. An 
acceptable compromise may be difficult to find and no systematic way 
exists to achieve it, even for this low order system. It is therefore 
obvious that more powerful methods are needed. 

The bandwidth of the system is determined by the selection of the 
weighting matrices A and B. Theoretically, state feedback controllers 
can have infinite bandwidth but in practical systems, the desired band- 
width will be limited by considerations such as saturation and noise 
susceptibility. It is of interest to compare the sensitivities of 
controllers with the same structure but with different bandwidths, in 
order to determine whether sensitivity considerations also contribute to 
the selection of the desired bandwidth. For this comparison, the weight 
ing matrices were changed and systems with two different bandwidths, 
one higher and one lower than the nominal, were calculated. Q and R 
were held constant, therefore keeping the estimator unchanged. The 
stability regions for these systems are shown in Table III-l. 

Table III-l 

FREQUENCY MARGIN AS A FUNCTION OF SYSTEM BANDWIDTH 


System 

Bandwidth 

, -I, 

(sec ) 

5u> 

l-i 

Qco 

u 

Low 

0.32 

-1.15 

+0.9 

Nominal 

0.8 

-1.0 

+0.8 

High 

1.8 

-0.65 

+0. 65 



An increase in sensitivity is observed as the. bandwidth is increased. 
Although this effect is not drastic, it is an additional factor that 
limits the system bandwidth that is achievable in practice. 

D. 6 Sensitivity Comparison of Different Controllers 

In this section the sensitivity of the state estimate feedback (SEF) 

controller is compared with that of two other controllers that may be 
used with this plant. They are: (a) a state feedback controller; 

(b) a classical network compensator. The first controller is only 
realizable if all the states are measured; which; in general; is not 
practical. It is used here only in order to emphasize the increase 
in sensitivity due to the use of SEF instead of state feedback. 

(1) State feedback . For this case, the characteristic equation 
is the determinant of (si - F - 5F + GC), viz., 

4 C 4 3 / C 3 \ 2 C 2 °1 

s + — s + I— + kls + — alts + — ak 

2 V 2 / 2 i 2 


. / 2 "2 ^1 

+ Is + — as + y a 

' 2 2 

The root locus as a function of &k is shown in Fig. II 1-8. 

From Fig. III-8 it is clear that no sensitivity problem exists for this 
case. The SEF controller for this system, while equivalent to the 
state feedback controller in many aspects when the system parameters 
are at their nominal values, is much more sensitive to parameter per- 
turbations . 

(2) Classical compensation . There is no unique compensation 
network for this system but any network must have the following char- 
acteristics: (i) provide lead at the first zero crossing of the 

system (point A) so as to have an adequate phase margin. (ii) provide 
an overall phase lag greater than 180° at co = in order to satisfy 

the stability requirements of Eq. (3.18). The simplest network which 
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FEEDBACK CONTROLLER 



has these characteristics is 


G (s) 
c 


(s + a) 
(s + to) 2 


(3.19) 


where a, to < 

In order to determine the values of the parameters c q , a, and b, 
a parameter optimization program was used in which these parameters 
were adjusted so as to minimize a cost function of the form of Eq. (2.3) 
with the same A and B matrices as were used for the design of the 
optimal controller for this system (Eq. 3.16), The result is: 


G (s) 
c 


3054 (s + 0.23) 
(s + 3. 2) 2 


(3.20) 


Both the nominal response and the sensitivity of this compensator are 
compared with those of the SEF controller. 

The closed loop eigenvalues of the system with the classical com- 
pensator are: -0,07 ± 4.85J; -0.81 ± lj ; -4.85; -0.3. Note the low 

champing of the first eigenvalue. 

For time response comparison, the responses to a step command 
and to a step disturbance were computed. These responses are shown in 
Figs. III-9 and III-10. The lower damping of the natural frequency 
in the classical compensator can be seen in the velocity and accelera- 
tion responses, especially in the disturbance response. For sensitivity 
comparison, the movement of the root at -0.07 + 4.85j as a function of 
Sk/k is shown in Fig. III. 11. All the other roots move towards stable 
zeros and therefore cause no instability. Comparing Fig. III-ll and 
Fig. III-4 , it can be seen that the classical compens'- tor (CC) is much less 
sensitive than the SEF compensator. The stability regions are 
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FIG. III-10 RESPONSE TO STEP DISTURBANCE OF STATE ESTIMATE FEEDBACK (SEF) CONTROLLER 
AND CLASSICAL COMPENSATOR. 






FIG. III-I1 SENSITIVITY ROOT LOCUS OF CLASSICAL COMPENSATOR 
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The classical compensator has very low sensitivity to increase in cu. 

This can be understood by considering its Bode plot, Fig. II 1-12. In 
this plot it can be seen that the phase angle reaches 180° only below 
w , whereas for the SEF controller, it reaches this value both below 
and above this frequency (see Fig. III-5). For the classical controller, 
the phase angle cp c is always less than 180°. 

E. CONCLUSIONS 

In this chapter the sensitivity to parameter perturbation of the 
SEF controller was analyzed and demonstrated by means of an example. 

In the specific example that was considered, the classical compensator 
is probably a better choice than the SEF controller since its lower 
sensitivity seems more important than its less acceptable time behavior. 
If, however, the full system and not just its low frequency approxima- 
tion is considered, it was shown by Bull [BU-1 ] that such a classical 
compensator will not provide adequate control unless the system bandwidth 
is lowered considerably. For the required bandwidth, SEF is hard to 
replace. 

The sensitivity problem for this full system is, however, just 
as severe as for the low frequency approximation. It is therefore 
clearly desirable to have a general method for the sensitivity reduction 
of SEF controllers. This method should operate in such a way that 
while the sensitivity is reduced, the nominal performance is not degraded 
unduly. 

In Chapter IV such a method is developed. 
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IV. 


A DESIGN METHOD FOR MINIMIZING THE 
SENSITIVITY TO PARAMETER VARIATIONS 




and (b) feedback control. If the control is implemented in an open 
loop fashion, 3u/d|l == 0, and the optimal control can be obtained in 
a straightforward way by solving the optimal regulator problem for the 
augmented system with the state vector x A = [x , cr ] , and the PI of 
Eq, (4,2). This, however, is hardly a realistic approach since, in 
general, feedback control will be required. 

In order to determine du/ci|_l for the feedback case the form of 
the control has to be stipulated. Some authors [CAS-1, DA-1, DO-1, BR-1 
stipulate 

u = Cj x + C^cr 

and solve the optimal control problem for the augmented system, 

and C >2 are not obtained from the solution of a Ricatti equation 
since the dynamic matrix of the augmented system contains C 
[SA-1 ] . 

In addition to some theoretical questions as to the optimality 
of this solution [SA-l], it is complicated to implement since the 
sensitivity vector cr has to be obtained from a model. The order of 
the system is thereby augmented. For single input system algorithms 
have been developed to augment the order by n only, even if there are 
several variable parameters in the system [WI-1], The examples given 
for the use of this method are of low order and assume that the state 
is available for feedback [CAB-I, DA-1, DO-1, LA-1, BR-1], No mention 
is made of the influence of using the estimate instead of the state. 
These examples seem to show a definite reduction in their trajectory 
sensitivity but it is not clear whether this is also true for higher 
order systems. In one case of a higher order system, no conclusive 
result was found [RY-2] . Hendricks and D’Angelo [HE-l] use the same 
augmented system but postulate a control of the form 

u st Cx 

and find C by parameter optimization. This method was applied to 
the sensitivity minimization of a space booster with good results. 
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Here too, however, all the states are assumed available and the effect 
of using state estimates instead of states is not considered. This 
and the following methods have the advantage over the previous methods 
that the controller does not become more complicated because of the 
sensitivity requirements and only the values of its parameters are 
modified. Rillings and Roy [RI-l] postulate the same control but use 
analog computer simulation to minimize the PI. In another variant of 
this method, Cassidy and Roy [CAS-l] force the control to have non- 
zero feedback gains only for the measured states. This is done by 
solving the inverse problem and defining the weighting matrices that 
give the desired feedback structure. From the given data, the compu- 
tation times for this method seem extensively lengthy. More recently, 
Stravroulakis and Sarachik [ST-l], using the same augmented system, 
derive iterative governing equations for output feedback or state 
estimate feedback gains for both the deterministic and the stochastic 
case. These equations are formulated for a single variable parameter 
and it is not clear how they can be extended to multiple parameters. 

The gains that are obtained by this method seem to be applicable to 
practical state estimate feedback systems but the governing equations 
are complicated and the computational labor involved in the iterative 
calculation maybe considerable. This problem is not discussed by the 
authors. Many of the references cited above and other papers that treat 
various aspects of the sensitivity problem were collected by Cruz in 
a book published recently [CR-1] . 

A different approach is used by Palsson and Whittaker [PA-1'1 . 

In this approach, the variable parameters are considered as components 

of a random vector with mean at the nominal value and known covariance. 

A performance index of the f orm . ... 

/ •. • ' / ■■■■ ■ -oo. -■ ' / 

J = ^ x T Axdt 

•; ; : , ••■'./. O ••• .-V •, ; : 

is minimized where the average is taken over the values of the variable 
parameters. The structure of the system is predetermined and the mini- 
mization is performed by varying a set of free parameters selected by 
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the designer. The method ns presented is applicable to single-input 
single-output (SISO) systems only, and is essentially designed to select 
the parameters of classical compensators. It seeks a balance between 
the performance at nominal and off-nominal conditions. The relative 
importance of these two conditions can be determined by the designer. 

The example presented shows an appreciable reduction in sensitivity 
without unduly affecting the nominal performance. This method is not 
suitable for the sensitivity reduction of SEF controllers for several 
reasons: (a) The natural free parameters for a SEF controller are the 

feedback and estimator gains. However, since the control is not weighted 
in the PI, not all the feedback gains can be left free since very high 
or even infinite gains would result in most cases. The selection of 
the gains that remain fixed is somewhat arbitrary and may lead to un- 
satisfactory results, (b) The restriction to SISO systems prevents 
the selection of the estimator gains as free parameters, if the estim- 
ator is a Kalman filter. In that case, at least two inputs are required; 
one being a disturbance, and the other a sensor noise (Eq. 2.36). If 
those two inputs cannot be used in the desensitization, the nominal 
estimator gains will tend to infinity and the estimator may lose its 
filtering properties. The total number of inputs must therefore be 
equal to the sum of the number of disturbances and outputs. The 
method as presented, however, cannot be extended to multivariable sys- 
tems. (c) The controller has to be transformed into its equivalent 
compensating network form or alternatively, the closed loop transfer 
function has to be found. Both these operations are cumbersome and 
present., numerical problems . 

In this thesis a sensitivty minimization method is developed that 
is based on the Palsson-Whittaker [PA-1] method but has none of the 
drawbacks described above. It is applicable to multivariable systems , 
both deterministic and stochastic, without restriction on the struc- 
ture of the system matrices. Although it was motivated by the need 
for sensitivity reduction of SEF controllers^ it is applicable to any 
system which can be represented in the form of E f. (2.1). The basic 
equations are similar to those given by Palsson-Whittaker but the met hod 
of solution is totally different. 



B< description of the sensitivity minimization method 

B. 1 Problem Statement 

Consider the system (repeated from 2.l) 


x = Fx + Gu + Tw 
y = Rx + v , 


(4.3) 


where the matrices F, G, I\ and H contain parameters the values of 
which are uncertain. If these parameters are considered as Gaussian 
random variables, a Gaussian random vector may be formed of which they 


are 


the components. This vector is specified by 


E(z) 


J n» 


(4.4) 


where the components of z n are the nominal values of the parameters, 
and by 


E[(z - z )(z ~ O 5 


(4.5) 


n 


a covariance 
be written as 


matrix which is assumed known. Equation (4.3) may then 


x = F(z)x + G(z)u + r(z)w 


H(z)x + v, 


(4.6) 


x(0) = 0. 


The control is defined as 


u = -Cx 


• (4.7) 


where the values of C may be left free or defined by functional re- 
lationships to other parameters of the system. For this definition of 
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u to be valid, Eq. (4.3) generally is required to describe an augmented 
system that includes the plant and the compensations. The matrices 
F, G, H, and T then have to be defined accordingly. For SEF controllers, 
the augmented system is given by Eq. (2.31). For this case, the control 
is given as 

u = [0, -C E ][x, x] = X -^ ' 


A quadratic PI for this system is 



lim 

t f -*» 


\ E(x T Ax + u T Bu)dt 


o 



(4.8) 


In this expression, the expected value is taken over the probability 
distributions of x and u that are derived from both the distributions 
of the random process w and of the random vector z. Note that w is 
the process noise of the augmented system. Since w and z are independ- 
ent, the expected value of a function of x and u is 


E[f(x,u) ] = ^ f [x(w,z), u(w,z)] p[x(w,z), u(w,z)]dxdy 


Xf u 


~ ^ E z p(w)dw = ^ E w p(z)dz , 


w 


where 


^ g(w,z) p(z)dz 


(4.9) 


is the expected value over the distribution of z and 


E 


w 


\ g(w, 


z) p(w)dw 


w 


is the expected value over the distribution of w. 

A free parameter vector q is defined by the designer. This 
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vector consists of n parameters of the system matrices that can be 

q 

varied by the design©!*# 

The problem is now stated as follows: Given the system of Eqs. (4.6) 

(4.7) in which the system matrices are functions of a variable parameter 
vector defined by Eqs. (4.4) and (4.5), determine the value of the free 
parameter vector q so that the PI of Eq. (4.8) is minimised. 


B.2 Method of Solution 

Using Eq. (4;7) in Eq. (4.8), the PI becomes 

% 

J = lim — — ( E[x (A + C BC)x] dt 


Since for any vector v, v T v = tr(vv T ), this equation can be transformed 

into tf 

X f T 

J = tr[( A + C T BC) lim — \ E(xx )dt] . (4.10) 

t J 


The state vector x may be written as 


x(w,z) = x n (w,z n ) + &x (w , 5z ) 


(4.11) 


where x is 
n 


the state vector obtained when the parameters have their 
nominal values, and 5x is the perturbation in the state vector due to 
a perturbation 5z in the variable parameter vector. 

Assuming small perturbations in z so that a first order expan- 
sion is satisfactory, substituting Eq. (4.7) into (4.6) and defining 


F = F - GC 
c 


the governing equations for and Sx become 


-79 




x = F (z )x + 

n c n n 

r(z )w 
n 


(4.12a) 


5x = F (z )5x + 
c n 

SF x + STw , 
c n 


(4.12b) 


x (0) = 0; 

n 

&x(0) = 0 ; 


(4.12c) 

where 

n z 

dF 




5F = £ 

C x^i 

c 

-s — Sz. 

CJz . X 

X 




n z 

sr - 2 

i-l 

dr c 

5z. . 

dz. 1 

X 



The expected value of Eq. (4.12b) is 





n z 3 >f 

n z 

dr 


E(5x) 

= F (z )E(5x) + Y\ — E(5z.x ) + 

c n i~l dz. 1 n i=l 

X 

— — E(5z.w) 
dz 1 ’ 

X 


5z 

is independent of x^ and 

w, and since E(Sz) = 0, 

we have 


E(6z. x ) = E(5z.w) = 

x n x 

0, i = 1 

t n z • 



From Eq. (4 (,12c) , 

E[5x(0>] = 0 

and therefore 

E[5x( t) ] =0. ■ (4.13) 

The expression E(xx^) in the PI of Eq. (4.10) can now be eval- 
uated in terms of x r and 5x 

E(xx T ) = E[(x n + Sx)(x n + 5 x) T ] 

(4 . 14) 

- E(x x T ) + E (x 5x T ) + E(6xx T ) + E(5xSx T ) . 
n n n n 
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To find the value of the second term of this expression, it can be 
written according to Eq. (4.8) and using Eq. (4.13), as 


E[x (w)Sx T (w,z) ] = E [E Cx (w)Sx T (w,z) ]} 

n w z n 

= E fx (w)E[5x T (w,z) ]} = o . 


w'- n' 


Similarly, 


We therefore get 


E[ & xx n ] = 0 , 


E(xx T ) = ■ E(x^x^) + E(6x5x T ) = X^ + SX , 


(4.15) 


where X is the covariance matrix of the nominal state. Sx can 
n 

thus be interpreted as the addition to the covariance due to the para- 
meter uncertainties. Substituting Eq. (4.15) into Eq. (4.10) yields 


J = tr[( 


t*£ 

A + C T BC) lim ( [X (t) + SX(t)]dt 

t -*» ■fcf ) n 


(4.16) 


= tr[( A + C T BC)(X n + SX) ] 


where X and 5X are the time averages over all time of X and SX 
n n 

respectively. 

For a stable system, X and SX tend to constant values X (*>) 

’ n n 

and SX (°° ) . Since the averaging in Eq. (4.16) is performed over a large 

time interval, it can be assumed that X -* X (°°) and SX -» SXC°), 

n n 

The PI therefore becomes 


J = tr{[A + C T BC][X («) + SX(»)]} 


n 


(4.17) 


" J 0 + J A ’ 
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where 


J q = tr[(A + C T BC) X n (oo>] 

is the nominal PI , and 

J = tr[(A + C T BC)5 X(o°)] 

A 

is the additional PI due to the parameter variations. The PI of Eq. 

(4.17) can be minimized computationally by a two-step sequence, (a) 
the matrices X n (») and SX(«0 are found for a given value of the free 
parameter vector, q. The PI that corresponds to this value is then 
determined. (b) the q vector is modified in a direction that de- 
creases J. This sequence is repeated until the decrease in J in 
one cycle is less than a predetermined value. 

The two parts of this sequence are independent and computer programs 
for each one can be developed separately. 

B . 3 The Governing Equations for X^ and 5X . 

For brevity, the subscripts of x and X^ will now be dropped. 

rp 

To find X use (d)/(dt)(xx ), viz., 

d , T N • T • T ( a -i o \ 

— (xx ) = xx + xx . (4.ia; 

dt 

Substituting Eq. (4.12a) into Eq. (4.18) yields 

T T 

(F x + rw)x + x(F x + rw) 
c c 

(4.19) 

T T TT T T 

F xx + rwx + xx F + xw r . 
c c 

Taking the expected value of both sides of this equation and using 
the definition 

e(xx T ) ; = .x 

Eq. (4.19) becomes : 
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d 

dt 


T 

(xx ) 



X 


(4.20) 


= F X + TE(wx T ) 4- XF T + E(xw T )r T . 
c c 


But , 


F c t 


F„(t-T) 


E(xw A ) = E{[e " x(0) + ^ e “ rw(T)d T ]w (t)} 

O 


and 


therefore 


e[w(t)w (t)] = QS(t - T) 
E[x(0)w(t)] = 0,for t ;> 0, 


E(xw T ) 


= E 


^ e F c( t "' r ^p w ( t)w T ( t)dT 


f F c (t-t) 

= 3 e c rQ&(t-T)dT 


= 2 r Q • 


Similarly, 


T 

E (wx ) 


1 T 

- or 

2 ^ 


Using Eqs. (4.21) and (4.22) in Eq. (4.20), the covariance is 


x = f x + xf t + rQr T . 

c c 


In the steady state, X = 0 and the final equation for X M is 


f x + x f t + rQP T = 0 . 

C » 00 c 


The governing equation for SX can be found in a similar way. 


(Sx&x T ) 

dt 


T T 

5x5x + Sxgx 
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(4.21) 

(4.22) 

(4.23) 

(4.24) 
Using 

(4.25a) 



and Eq. (4.12b) yields 


H T 1 TT TT TT 

— ( BxBx ) = BxBx F + Bxx 8F + Bxw BP 

dt c c 

+ F Bx8x T + BF xfix^ + BrwBx^ . 
c c 


(4.25b) 


To determine the expected value of this equation, E(8xw ) has to be 
evaluated. 


E( Bxw^) 


F c t 
E I f e 8x 


■ s 


( c 

t P C (t-T) 


o + ^ e F °^ T ^[BF x(t) + BPw(t) ]5t}w T ( t) 


BF E[x(t)w (t) ]dT 


(4.26) 


ft F_(t-t) T 

+ V e BPE[w(t)w (t) ]dT 


o 


= f 5PQ . 


The first integral vanishes because 


E[x(T)w T (t) ] £ 0 only for T = t . 

Using the definitions of X and 8X from Eq. (4.15), the expected value 
of Eq. (4.25b) can now be found. 


SX = BXF T + F BX T + E[8rQ8P T ] 


c c 


+ E[8xx T 5F T ] + E [ BF xSx T ] . 
c c 


(4.27a) 


The last term on the right hand side of this equation may be written as 


E[SF„xdx T ] = e z [8F c E w (xBx T ) ] , 
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since 5F is independent of w. Defining 
c 

Y = E w (x8x T ) , (4.27b) 

the 'Last term of Eq. (4,27a) is 

E[SF x5x T ] = E [5F Y] . (4.28) 

u c z c 

As t -» 5X tends to the steady state value, SX^. 

The final equation for SX^ is obtained by substituting Eq. 
(4.28) into Eq. (4.27a) and putting SX^ = 0. Thus 

SX F T 4- F BX + E [eTQST T + Y^5F T + 5F Y ] = 0 (4.29) 

OO Q C°5 z 00 C 

To get the governing equation for Y^ the identity 


“r (x5x T ) = x5x T + X&X T 
at 


m T T T T T T T 

= F x5x + TwSx + x&x F + xx SF + xw &r 
C c ° 


is used. The expected value of this equation over the distribution of 
w can now be found using Eq. (4.22) and Eq, (4.26) , and the definition 
of Eq. (4.27b). 

Y = F Y + YF T + rQ&r T + X5F^ . (4.31) 

The steady state value of Y as t -» » is given by 

F Y + Y F T + rQ5r T + X BF^ = 0 . (4.32) 

C eo 00 c oo c 
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The governing equations for and SX^ that have been derived 

in this section are rewritten 


F X + x/ = 
q oo oo c 

-ror T 

(4. 33a) 

T 

F Y + Y F 
c 03 00 c 

-rw T - x m sf^ 

(4. 33b) 

F 5X + SX F T 

q oo c» c 

= -E z [BrQSr T + VV + SF„YJ. 

(4 j 33c) 


These equations all have the same well known form [e.g., GA-l] 

T 

AX + XA = 3 , 


where A and B are given. They must be solved in the order in which 
they are written because the right hand sides of the second and third 
equations contain expressions found in the solution of the previous 
equations. For their actual solution, these equations must be written 
in a slightly different form. Equation (4.33b) is rewritten as 


F Y + Y F 
c 00 07 c 


- sfrsr 1 

i =1 \ 0z. 


3f 


- x “sr r* 
1 


(4.34) 


y is therefore obtained in the form of a sum: = £ Y,5z, where 

CO 

Y is the solution of 
i 


t ar BF c 

F Y. + Y.F T = - FQ - X 5—- . 

ci i c o z_^ m ® Z i 


(4.35) 


This equation is solved n times with different right hand sides, 

z 

The right hand side of Eq. (4.33c) can be written as 



E [srQ&r T + SF Y + Y^8F T ] 

z C 03 00 c 


= E 


n z n z \T Sf /Sf 

Lfi £ 51 ' “(fc)* s;v , 


T_ 

/Sf \ 

lv 


Sz i 8z j (4.36) 


" S. [* /ar V + _ + / ap c _ V 

" £ £ ST Q lc 5 T J + ST Y J + Y j) 

j=l x=l |_ i ' 3 ' 1 x ' J 


v . . 
1 3 


where V is the parameter covariance matrix (see Eq. 4.5). 

Equation (4.36) can be written in a more compact form although 


this form is not used 

in the 

storage requirements. 


n z n z r> 

E Ef 

/Sr \ T 

Q s~) 

j=i i=iL az i 

va V 

rsr i r ' T 

1 Sr T 

= [^<*][ v ®i 

npj S z 


f Sf \ t] 




(4.37) 


1 X" T + X’ ( V ® 1 )Y + Or (V ® i)y]’ 
dz Oz n |_S Z n J 


where I is the n x n unit matrix, 
n 


Sr 


Q = 


Sr ! Sr „ 

7“ Q ! -n Q 

Sz i i 3z 2 


Sr 


Q 


iz 


i [n x(n X n )]; 

P z 


Sr 

S z 


Sr 


L i 


Sr 


[n x(n X n ) ] ; 

p Z 
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Bf 

c 


C>F 

c 

I 1 

1 1 

t • • • 1 

Sf 

c 

dz 


>1 

1 1 
1 1 
1 1 

Sz n 

U .Z J 



F- — 

Y 

1 



Y 

= 

• 

• 

• ■ 





Y r 

r *Z 




[n X (n X n ) ]; 

z 


■[(nXn)xn] j 


and where 


A ® B 


the Kronecker product of A and B, This product is defined for 
A[n X n] and B[m X m] as a (n X ra) X (n X m) matrix having as 
its ijth m x m block: 


(A <g> B) = a . . B. 
1 3 


There are n X n such blocks. 

Equation (4.33) can now be rewritten in the form 


T T 

F X + X F = -ror 

C 00 oo c ^ 


F Y + Y.F 
c x i c 


F 5X + 8XF 
c 00 00 c 


-I Q * - ^ i 
3z. 

i l 


-MM© 

- — (v ® I )y - (v ® I )y 1 
C> z n bz n J 


(4. 


(4. 


(4. 


38a) 


38b) 


38c) 
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From Equation (4.33c) and (4.36) it can be seen that 5X is 
proportional to the parameter covariance matrix V, i.e., if each 
element of V is multiplied by a scalar coefficient e, so is each 
element of SX. From Eq. (4.17) it is clear that the additional per- 
formance index J A is then multiplied by the same coefficient. 

In the problem statement the parameter covariance matrix V is 
assumed known (see Eq. 4.5). Its elements are a measure of the un- 
certainties in the parameters. In reality, however, these uncertain- 
ties are ill defined and the matrix V is considered mainly as a de- 
sign tool. It can be written as 


V = 6V 0 

where c is a scalar weighting coefficient. Equation (4,17) then 
becomes 


= J 0 + J A 


J 0 + eJ A0 


(4.39) 


The coefficient e is selected by the designer for the relative weight- 
ing of nominal performance vs sensitivity. The relative magnitude of 
the elements of V is selected according to the importance of the 
sensitivity reduction for the respective parameters. 

Equations (4.38) remain valid if the system is forced by inputs 
other than white noise. If the disturbance is colored, it can be 
modeled by means of a shaping filter forced by white noise. An aug- 
mented system can then be formed consisting of the states of the system 
and of the shaping filter. This augmented system is excited by white 
noise and Eqs. (4.38) is therefore valid for it. 

For deterministic systems which are required to recover from non- 
zero initial conditions, the white noise vector w is replaced by an 
impulsive input vector w. at t = 0: w^ = w^S(t), The initial con- 

ditions are then represented by the equivalent initial impulses . 
Equations (4.38) can be used unchanged if the following terms are 
redefined: 
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Q 


T 

= w w 
o o 


^ xx^dt 

Y oo = \ X5 * T dt 

■o 

CO 

8X = ( SxSx^dt 

00 ) 

J o 


(4.40) 


With these redefined terms the deterministic PI is given by Eq. (4.17). 

If the system is to be optimized for some deterministic input 
other than an initial impulse, this can be handled by state augmentations 


B.4 Description of the Computer Program 

The computer program, PAROPT, for the minimization of the PI 
of Eq. (4.7) is described in detail in Appendix B. Only its main 
features are described in this section. It consists of two main 
parts: (1) A search subprogram; (2) A subprogram for the solution 

of Eq. (4.33) and determination of the PI from Eq. (4.17). 

(1) Search subprogram . This subprogram is part of the program 
library of the Computer Science Department at Stanford 
University. It was developed by Gill, Murray and Pitfield 
[GI-1 ] . It iteratively seeks the minimum value of a scalar 
function F(q), where q is a vector of dimension n 
by modifying the components of q. The value of F(q) for 
a given q is an input to the subprogram. For the current 
problem, it is the performance index J of Eq. (4.17) . 

The mode of operation of this program is as follows. For 
given initial values of q and J, the approximate direction 
of the gradient of J with respect to the vector q is 
found by perturbing the components of q one at a time and 
determining J for the perturbed vector. The components 
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of the approximate gradient vector are proportional to 


The direction of the conjugate gradient [BRY-2 for definition! 
is then determined and a linear search is performed along 
this direction. In general, three to four evaluations of 
J are required to determine the minimum of J along a direc- 
tion. At the minimum, the gradient is again found and a 
new linear search direction is determined. 

One iteration is defined as a gradient determination and 
linear search. This iterative procedure is continued until 
termination criteria are satisfied. (See Appendix B for 
more details.) 

(2) Subprogram for the evaluation of J . It is obvious from the 
description of the search subprogram that the value of J 
must be computed a large number of times. In one iteration, 
n evaluations are required for the gradient determination 

q 

and 3 to 4 for the linear search. In an average program, 

8 to 10 iterations may be expected and therefore it may be 
required to compute J 100 to 200 times. For the program 
to be of any practical usefulness, a very efficient method 
for this computation must therefore be developed. 

The principal part of this evaluation is the solution of 
Eqs (4.38). The same equation with different right hand 
sides is to be solved 2 + n^ times. It is to be noted, 
however, that Eq. (4.38a) and (4.38c) have symmetric right 
hand sides and therefore symmetric solutions (Lyapunov 
equations) , whereas the right hand side of Eq. (4,38b) is not 
symmetric. Several methods are available for the solution of 
these equations. These methods are compared by Hagander [HA-1] 
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and Pace and Barnett [PA-l], The recommendation of Hagander 
is to use a direct method (described below) for systems 
of order smaller than 6 or 7, and other methods above this 
order. The reason for the recommendation for a different 
method for large systems, despite the better precision of 
the direct method, is that the computation time of the 
direct method increases as n®/3 and becomes prohibitively 
large for large systems. Pace and Barnett recommend the 
direct method up to order 10 approximately. The direct 
method consists of transforming an equation of the form 

AZ + ZA T = B (4.41) 

into the form 

CLz = -{3 , (4.42) 

2 

where z(JL X 1) is the n vector of the coefficients of 

2 

Z, and |3 is the n vector of the coefficient of B. 
a is obtained from A by [BE-1] 

a = A ® I + I ® A 
n n 

For symmetric Z, the dimension of z need only be 
H = (n + l)n. An algorithm to obtain CL for this case 
is given by Bingoulac [BI — 1 ] . 

For antisymmetric Z, the dimension of z is l = |(n-l)n. 

A transformation algorithm for this case was developed as 
part of this program. 

The method commonly used for the solution of linear equa- 
tions of the form of (4.42) is Gaussian elimination [F0-1 ]. 

It consists of two steps: (a) forward elimination CC /3 oper 

ations), and (b) back substitution J 2 operation). By 

comparing the number of operations, it is obvious that to 
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solve a general equation of type (4.41) by the direct 
method, it is advantageous to decompose the right hand 
side into symmetric and antisymmetric parts and to obtain 
two separate solutions, the sum of which is the required 
solution. This is important for the solution of Eq. 

(4.38b). In the solution of (4.42), forward elimination is 
only required if a is changed. If only P is changed, 
the much less costly back-substitution is required. 

In our case, for one evaluation of the PI, only one forward 
elimination of order Kn+l)n is required (for X), and one 
of order £(n-l)n (for the antisymmetric part of . This 
considerably reduces the average computation time for the 
direct method and makes it attractive for much larger sys- 
tems. Moreover, since the points at which the cost is 
computed for the gradient determination correspond to only 
slightly perturbed q vectors, the changes in X and SX 
that stem from these perturbations can be obtained by ex- 
panding Eqs . (4.38) about the nominal point and neglecting 

second order terms. Equation (4.38a) at the perturbed 
point is 


r x« r I 


where the subscript 1 refers to this point. Expanding 
Eq. (4.43) yields 

: (F +DF ) (X+DX) + (X+DX)(F c +DF c ) T = 

where the prefix D indicates the change between the nom- 
inal and perturbed points. Neglecting the product (DI< C )(DX), 
Eq. (4.43) becomes 


F 

c 





T T 

xf , - r qt 

cl 1 1 


(4.44a) 


The coefficients of the left hand side of this equation are 
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the same as for Eq. (4.38a) at the nominal point. Similarly 
for DY and D(SX), 


F DY + DYF T = 
c c 


F D(SX) + D(SX) 

C 


V + YF Il ' r i Q6r I ' Vol 


= -F BX - BXF 

cl cl 


-E ( 5F Y + Y.5F T . + 5T QBT^) 
z cl 1 1 cl 11 


(4,44b) 


(4.44c) 


This approximation was compared to the exact method of 
evaluating the perturbed PI and no significant difference 
was observed. 

The evaluation of the PI and its gradient, therefore, can 
be done with only one forward elimination of order ^(n + l)n 
and one of order ^~(n — l)n. This method was implemented in 
PAROPT, which has been applied to the sensitivity reduction 
of several systems. Results of its application are described 
in Section IV--C, 

The largest system to which it was applied is of 12th order 
with 2 variable and 20 free parameters. 

One iteration for this system (gradient evaluation plus 
linear search) required about 40 sec on an IBM 360/67 
computer. This method therefore seems practical for even 
fairly large systems. The computation time increases as the 
fourth power (approximately) of the system order and is almost 
independent of the number of free and variable parameters. 

C. APPLICATIONS 

c.l Introduction 

In this Section the application of the sensitivity reduction 
program ( SRP) to two systems will be described. The systems are: 

( a ) low frequency approximation of the Stanford Relativity Satellite (SRS) 
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as described in Section III-D. (b) full Stanford Relativity Satellite 
(per III-D). 

The design considerations for the nominal optimal controllers 
for these systems will not be given since they are described in detail 
by Bull [BU-l] for the full system and remain valid for the low fre- 
quency approximation. Since the purpose of this section is to investi- 
gate the operation of the sensitivity reduction program (SRP) , several 
designs with different program parameters were made and the results 
are presented in considerable detail, particularly for the first exam- 
ple. The criteria that are used for the comparison of different de- 
signs for the same system are: 

• Sensitivity criterion: the range of variation of the variable 

parameter for which the system is still stable (high range = 
low sensitivity); 

• Nominal performance criteria : output and control rms values 

and the square root of the nominal PI, which is a weighted rms 
value. 

Some points have to be kept in mind when using these criteria: (a) 

The stability range of the variable parameters should not be construed 
as defining the actual permitted range of variation of these parameters. 
In general, the performance will become unacceptable for variations that 
are considerably less than those that cause instability. (b) The rms 
values of the outputs and the controls depend on the assumed covariance 
matrices of the process and measurement noises. These covariance 
matrices are generally not well known and in some cases, they are 
artificially determined in order to get acceptably damped roots of the 
estimator. Small differences (less than a factor of 2) in the rms 
values of different designs cannot, therefore, be considered as signif- 
icant. The criteria that are used in this section are therefore un- 
refined but can still give a valid comparison between different designs. 

More precise criteria are difficult to define, in general, al- 
though for specif ic cases they may exist. In some cases the time response 
envelope to a specific input may be restricted, or limits may be posed 
on the phase and gain margins. It is important to verify, whenever 
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such criteria are used, that they reflect actual system requirements and 
do not pose artificial restraints on the design. 

The initial point in the application of the SRP must be a stable 
one for Eq. (4.33) to be valid. In the examples of this Section, the 
nominal optimal point that was found without considering sensitivity 
was selected as the initial point. This is not required but it is an 
easy point to calculate since the weighting and covariance matrices are 
also required for the SRP. The systems that are used in the SRP are the 
augmented systems given by 



(4.45) 


where F is the actual plant, F is the assumed plant in estimator. 

At the initial point F = F. The weighting matrices are the same 
as those that were used for the nominal design. 

For the application of the program, the parameter covariance 
matrix and the sensitivity coefficient 0 have to be determined and 
the free parameters have to be defined. The parameters that are con- 
trolled by the designer and some or all of which can be used as free 
variables are: C, the feedback gains; K, the estimator gains; and 

F, the representations of the variable parameters in the estimator. 

The last item may require some clarification. If there are no variable 
parameters, the estimator parameters will obviously be selected to 
be the same as the plant parameters. If, however, some plant parameters 
are variable, it has been found that the sensitivity can be decreased, 
in some cases, if their representations in the estimator differ from 
their nominal values. It is therefore desirable to include these repre- 
sentations among the free parameters. The actual selection of the pro- 
gram parameters will be discussed for each example separately. 




C-2 Example 1 : Low Frequency Approximation of the Stanford 

Relativity Satellite (SRS) 

This system was described in. Chapter III-D. The augmented system 
is an 8th order system lor which the initial dynamic and state weight 
matrices are shown in Pig. IV- 1. 

(1) P rogram parameters . The only variable parameter for this 

case is the spring stiffness, k. The parameter covariance 
matrix is therefore a scalar and only its product with the 
sensitivity coefficient £ is important. Various values 
of this product have been used in the different designs 
as described below. Only the estimator gains and the vari- 
able parameter representations in the estimator were used 
as free parameters since it was found in preliminary runs 
that the feedback gains do hot vary appreciably. The 
application of the SRP (sensitivity reduction program) in 
this case is therefore a redesign of the estimator. 

Four designs were executed with the design parameters 
varied as described below. The weighting coefficients 
6 and V (a scalar in this case) define the assumed rms 
value of the spring constant k 

CTj, = ^ev' (4.46) 

As explained in IV-B-3, this rms value does not represent 
an actual expected uncertainty in the variable parameter 
but is used as a design tool for the relative weighting of 
the nominal and additional performance indices. The designs 
are: 


Design No. 1: 

<T k /k = 0.28 (J a /J q = 1 

at the initial point) 

Design No. 2: 

<T k /k = 0.9 (J a /J q = 10 

at the initial point) 

Design No. 3: 

cr k /k = 0.9, R = 0 


Design No, 4: 

cr /k = 0.9, Q = 0. 

k 
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FIG. IV- 1 DYNAMIC AND WEIGHTING MATRICES FOR THE REDUCED STANFORD 
RELATIVITY SATELLITE 



In Designs 1 and 2 , the estimator is a Kalman filter and 
the SRP seeks a balance between the nominally optimal 
estimator parameters and those required for minimum sensi- 
tivity. In Design 3 it is assumed that there is no measure- 
ment noise. The minimum nominal PI is obtained for K -» “, 
since for this gain, the estimate error covariance P -» 0 
and the system behaves as if state feedback instead of state 
estimate feedback were used. The estimator gains are there- 
fore determined by sensitivity considerations only. Simi- 
larly, in Design 4, where no process noise is assumed, the 
minimum nominal PI is obtained for K = 0. Here, too, the 
estimator gains are determined by sensitivity considerations 
only. 

(2) Results . The values of the estimator parameters for the 
different designs are given in Table IV- 1. These include 
the estimator gains and the value of the spring constant 

Table IV-1 

ESTIMATOR PARAMETERS OF THE 
REDUCED STANFORD RELATIVITY SATELLITE (SRS) DESIGNS 



Nominal 

Design 
No. 1 

Design 
No. 2 

Design 
No. 3 

Design 
No. 4 

k l 

18.7 

4.46 

31.1 

26. 3 

45.6 

k 2 

175.4 

166.2 

162.9 

163.7 

179.5 

k 3 

13.9 

19.8 

107.8 

108.7 

134.6 

k 4 

-30.5 

-111.6 

-50,5 

-49.6 

-34.4 

Assumed 
Value of 

k <“o> 

25 .0 

26.4 

22.3 

,7 " . 

■ . 

22.4 

31.6 
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used in the estimator (k). The nominal performance criteria 
are given in Table IV-2. These are the average, output and 
control ruts values that result from the process and measure- 
ment noises given in Chapter III-C. 

The eigenvalues are shown in Table IV-3 and in Fig. IV-2. 

The sensitivity root loci for the sensitive eigenvalue are 
given in Fig. IV-3. The frequency margin for the nominal 
design (from Fig. III-7) and for Design 2 are shown in 
Fig. IV-4. The sensitivity properties are compared in 
Table IV-3. 

Table IV-2 

NOMINAL PERFORMANCE CRITERIA OF THE REDUCED SRS DESIGNS 

(scaled to initial value) 



Nominal 

Design 
No. 1 

Design 
No. 2 

Design 
No. 3 

Design 
No. 4 

Nominal PI 
(J 

0 01 

1.0 

1.18 

1.59 

1.57 

1.32 

Weighted rms 
' o' Ol 

1.0 

1.09 

1.26 

1.25 

1.15 

Output rms* 

W 

1.0 

■ 

1.14 

1.72 

1.70 

1.74 

Control rms* 

(d /a ) 
u ul 

1.0 

1.12 

1.11 

1.15 

1.11 


* Note: The output and control rms values for all the designs were 

found using the same process and measurement noises 
(given in Sec. III-D-2) 




Estimator i Controller 


Table IV- 3 


EIGENVALUES OF THE REDUCED SRS DESIGNS 


Nominal 

Design No. 1 

Design No. 2 

Design No. 3 

— — *1 

Design No. 4 

-0.36 ± 5.01J 

-0.38 ± 4 , 98j 

-0.38 ± 5 . 05 j 

-0,38 ± 5. 05 j 

-0.32 ± 4 . 9j 

-0.41 ± 0,41 j 

-0.41 ± 0.41 j 

-0.41 ± 0.41 j 

-0.41 ± 0.41 j 

-0.41 ± 0,4 j 

-1.02 i 5. 11 j 

-1.41 ± 3. 17j 

-1.16 ± 9.19J 

-1.21 i 9.95J 

1 

o 

It 

I-* 

o 

w 

-8.33 ± 8. 35 j 

-0.8 ± 13. 2j 

-27.6 

-22.6 

-43.2 

• 


-1.14 

-1.21 

-1.08 


Jo (see l ) 



FIG. IV- 2 EIGENVALUES OF NOMINAL AND DESENSITIZED 
DESIGN OF REDUCED STANFORD RELATIVITY 
SATELLITE. \ 
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FIG. IV- 3 SENSITIVITY ROOT LOCI OF REDUCED RELATIVITY 
SATELLITE . COMPARISON OF NOMINAL AND DE- 
SENSITIZED DESIGNS. ' 
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FIG. IV-4 FREQUENCY MARGIN COMPARISON OF NOMINAL AND DESENSITIZED 


DESIGNS. 
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Table IV-4 


SENSITIVITY PROPERTIES- -STABILITY RANGE OF THE REDUCED SRS DESIGNS 



Nominal 

Design 
No. 1 

Design 
No. 2 

Design 
No. 3 

Design 
No. 4 

Range of 
Stability 

+ Ak/k 

+0.26 

+ 1.0 

+ 1.27 

+ 1.6 

+ 1.92 

-Ak/k 

-0.27 

-0.52 

-0.8 

-0.84 

-0.8 

Total 

Ak/k 

0.53 

' 

1.53 

1.52 

2.07 

2.44 


(3) Evaluation of the results . (i) As expected, increasing 
e decreases the sensitivity and increases the output and 
control noise. In Design 1, the stability range is in- 
creased by a factor of approximately 2.5 with an output 
rms increase of 14% only. If the stability range is in- 
creased by a factor of « 3.5 (Design 2), the output noise 
increases by 72%. This design, therefore, constitutes 
a considerable departure from the nominal optimum. Due 
to the somewhat artificial nature of the nominal optimum 
(see iv-C-1), this difference cannot be considered as very 
significant. A better balance between the increase in 
output and control noise can probably be achieved by trial 
and error. It is important to note that the differ- 
ence in stability range between Designs 1 and 2 does 
not fully account for the improvement of Design 2 over 1. 

In Fig. IV-3, it can be seen that in Design 2, the root 
locus does not approach the imaginary axis appreciably 
for a range of variations of k that is close to the 
stability range. Therefore, the practically acceptable 
range of variations (e.g. , | Re \\ > 0.2) is much larger 
for this Design than it is for Design 1. 

(ii) The lowered sensitivity of Design 2 can also be ob- 
served from its frequency margin (Fig. IV-4) which is much 
larger than that of the nominal system. 
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(iii) The estimator eigenvalues of Designs 3 and 4, 
particularly 3, are close to those of Design 2. This 
indicates that the estimator gains for Design 2 are deter- 
mined largely by sensitivity considerations, 

(iv) The two extreme cases represented by Designs 3 and 4 
indicate that for this system the estimator gains may be 
determined by sensitivity considerations solely, and still 
give adequate nominal performance. 

(v) Even for this simple system, it would have been 
difficult to determine the minimum sensitivity eigenvalue 
locations without some general method such as the one 
that was used. The removal of the eigenvalue at -1,02 

+ 5 . 11 j from the vicinity of the plant eigenvalue at 
-0.36 + 5 . 01 j looks plausible; but the movement of 
the eigenvalue at -8.33 + 8.35j is difficult to justify 
intuitively. 


C-3 Example 2: Full Stanford Relativity Satellite (SRS) 


(1) Program parameters . The augmented system is a 12th order 
system with two variable parameters--the stiffnesses of 
the two springs. The variations of these stiffnesses 
are unrelated and therefore the parameter covariance matrix 
has only diagonal elements. The design parameters that 


have to be determined are e, 


Different values of these parameters are used for the 
different designs (Table IV-5). The feedback gains, 
estimator gains, and assumed values of the spring constants 
are used as free parameters in all cases (a total of 20 
free' parameters) . ..x 
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FIG. IV- 5 DYNAMIC AND WEIGHTING MATRICES OE THE STANFORD RELATIVITY SATELLITE 



The initial dynamic and weighting matrices were shown 
in Fig. IV- 5. 

Three controllers were designed. They are described in 
Table IV-5. 


Table IV-5 

DESIGNS FOR THE STANFORD RELATIVITY SATELLITE 


Design No. 1 


0 /k 

V a 

1 

0.55 

0.18 

2 

2.5 

0.81 

3 

5.0 

0.81 






The motivation for these various designs is evident from 
the Table. As explained earlier, the parameter rms values 
are used as a design tool only for achieving the required 
sensitivity. 

(2) Results . The final values of the free parameters for the 
different designs are given in Table IV-6. The eigenvalues 
are shown in Table IV-7. The nominal performance criteria 
are compared in Table IV-8, The process and measurement 
noises used for the determination of these criteria are 
those used in Bull [BU-1], 

The sensitivity properties are compared in Table IV-9. 
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Table IV-6 


FREE PARAMETERS (in MKS units) OF THE FULL SRS DESIGNS 






Design 

Design 

Design 




Nominal 

No. 1 

No. 2 

No. 3 



C 11 

18391.0 


18391.0 

18391 . 0 


05 

c, „ 

632.8 


631.1 

631.2 



12 





| g 

S Ph 

C 13 

1881.8 

3 

a 

3 

1880.7 

1880.5 

CQ 

S5 H 
H O 

C 14 

126.1 

rH 

0 

o 

127.6 

127.8 

!5 

O 

c„ _ 

2285.2 


2285.4 

2285.5 

w 

<q 


15 

CQ 



O 


c „ 

-2001.3 

c 

-2001.6 

-2001.5 

w 

o 


1 6 

- g 




C 

101.4 

is 

99.7 

99.7 

< 

CQ 


21 




Q 

W 

w 

05 

W 

C 22 

274.8 

cn 

03 

254.9 

256.1 


05 13 

W 

Eh O 

C 23 

115.6 


70.9 

73.1 


B 05 
O Eh 
g 

C 24 

. • . 

279.8 

B 

CQ 

CA 

310.5 

310.8 


O 

O 

C 25 

275.6 


278.1 

277.9 



c 

372.8 


372.9 

371.3 



26 






k l 

24.1 

23.0 

24.2 

43.2 



k 

292,0 

297. 7 

287.2 

283.4 


05 

2 





O 

H M 

k 3 

-2.2 

-10.3 

-10.9 

-19. 1 


<5 55 

S M 
M <5 

k 4 

-124.0 

-205.4 

-219.0 

-224.0 


Eh O 
CQ 

w 

k 5 

14.9 

4.8 

-3.6 

-10.9 



k 6 

-26.0 

-41.5 

-41.5 

-53.6 

Plant 

Parameter 

mm 

625.0 

736.0 

737 . 0 

736.0 

Representation 

mSki 

25.0 

25.8 

33.2 

27.7 

In Estimator 

wm 
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Table IV- 7 


EIGENVALUES OF THE FULL STANFORD RELATIVITY SATELLITE DESIGNS 


NOMINAL 


DESIGN 1 


REAL 

-49,07)45521 
-49.07045521 
-0.39369490 
-0.39369490 
-0. 3051*4360 
—0 .385 14360 
-2.42507047 
—2. 4250 7047 
-8.62282277 
-8.62282277 
-1,05210676 
-1.05210676 


I MAG . 

48.32029403 
-48.32029403 
4.98784394 
-4.98784394 
0.47217532 
-C. 47217532 
2 5 *50344618 
-25.58344618 
7.12131077 
-7 « 1213 1077 
5. 11812296 
-5.11812296 


REAL 

— 4 S . 61510571 
-48. 61518571 
-4 . 09G83984 
-4*09983984 
-6.81374959 
-6.81374959 
-0,3436-9446 
-0.34369446 
-1.17993306 
-1.17993306 
-0.34520706 
-0.34520706 


DESIGN 2 


1 MAG . • 
49.93853088 
-49,93853083 
24,00645133 
-24.0064513-3 
7.84625316 
-7.84625316 
4.99512259 
-4.99512259 
4.158,78 479 
-4 .15878479 
0.40454027 
-0.40484027 


DESIGN 3 


REAL 

-48.10299920 
-48 • 10299920 
-3 *79056908 
-3^79056988 
-7 • 244 10743 
-7.24410743 
-C. 4 50 64515 
-0.45084515 
-1.61314669 
-1.6131466? 
-0.38407223 
-0. 38407223 


IMAG. 

49,32370542 

-49.32370542 

24.25442200 

-24.25442200 

4.30428075 

-4.30428075 

4.86150763 

-4.816150763 

3.48569047 

-3.48569047 

0.17419996. 

-0.17419996 


REAL 

-48.21998580 
-40.21998500 
-35.00701057 
-2.5406401C 
-2.54864010 
-0.41573592 
-0.41573592 
-0. 80529697 
-0.80529697 
-0.38388131 
-0.38380131 
-1.6109902.2 


IMAG. 

49.99571576 

-49.99571976 

0.0 

21.88676425 
-21.08676425 
5.03969535 
-5.03969535 
3.78422822 
-3.78422822 
0. 10814120 
-0. 18814120 

o.n 
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RANGE OF STABILITY 


Table IV-8 


NOMINAL PERFORMANCE CRITERIA OF THE FULL SRS DESIGNS 
(scaled to initial value) 



— 

Nominal 

Design 
No. 1 

Design 
No. 2 

Design 
No. 3 

Nominal PI 
[J /J .] 

u o Ol 

1.0 

1.02 

1.15 

1.22 

Weighted rms 

1.0 

1.01 

1.07 

s 

Output rms 

AA 

1.0 

1.27 

■ 

1.88 

■ 

Inner Control 
rms 

t>ul /0 uli ] 

1.0 

' ‘ 

1.28 

1.49 

1.44 

Outer Control 
rms 

KAsi 3 

— 

"T 

1.0 

0.99 

1.02 

1.02 

. . 


Table IV-9 

SENSITIVITY PROPERTIES OF THE FULL SRS DESIGNS 




— ; — 


Design 

Design 

Design 




Nominal 

No. 1 

No. 2 

No. 3 



+ '^ 7 / k 7 

+0, 71 

+0.97 

+1.3 

+ 3.1 

>- 

o 

T* 

♦—4 

a 

a* 

to 

fa 

- “A 

Total 

-0.12 

-0.22 

-0.24 

-0.23 

3 

m ■■ 

a 

< 

H 

M 

H 

CO 

A 

0.83 

1.19 

1.53 

1 ■ ; 

3.31 

o 

w 

p 

O 

Z 

+ “<A 

+0.32 

+0.8 

+1.4 

■ V- ■■ 

+1.2 

< 

a: 

to 

- ^ A 

i 

o 

w 

to 

-0.4 

-0.44 

-0. 38 


SOFT 

Total 



: • ' 




“A 

0.64 

1.2 

1,84 

1,58 



V J A 

1.0 

0.76 

0.7 

0.65 


nom 


; - ■ .• : . - 
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(3) Evaluation of the results . (i) The increase in the output 
and control rms for this example is much larger than it 
was for the Example 1 (compare Tables IV-2 and IV-8). 

This increase, however, is not reflected in the nominal 
PI (Table IV-8). This PI is a weighted average of the 
state and control rms values in which other states, beside 
the output, are weighted. If the rms of some of these 
states is reduced by the SRP, this reduction can balance 
the increase in the output rms. A better balance between 
sensitivity and nominal performance can probably be obtained 
if only the output and the controls are weighted in the PI. 

(ii) The range of stability is not symmetrical about the 
nominal value especially for the stiff springs. All the 
designs are more sensitive to the decrease of the spring 
constant of this spring than to its increase. In Table 
IV-6 it can be seen that the assumed value of this con- 
stant in the estimator is higher than the nominal. It 
seems, however, that for a more symmetric range, it should 
be lower. The stability range of a modified design, Design 
3A, was therefore found. This design has the estimator gains 
of Design 3 but the assumed value of the stiff spring con- 
stant in the estimator is 


2 ~2 

(The nominal value is w — 625 sec ) . The performance 

of Design 3A is shown in Table IV-10. Comparing this 
Design with Design 3, it can be seen that although its 
total stability range is lower, it still may be preferable 
since it is more symmetric. There is no obvious way to 
introduce such considerations into the program but they 
may lead to modifications of the program results by the 
designer. The nominal performance criteria are compar- 
able for the two Designs. 
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Table IV-10 


PERFORMANCE OF DESIGN No. 3A 




= 500 sec 

- 1 ) 

Range of 
Stability 

+Ak /k 
7 7 

+ 1.4 

For Stiff 
Spring 

-£k /k 
7 7 

-0.48 


Total 



Ak /k 

7 7 

1.88 



Output rrns 

E V a 9i] 

Inner Control rms 

r cr /a ] 

L ul ; uli J 

Outer Control rms 

C WW 



(iii) From Table IV-6, it can be seen that the feedback 
gains are changed very little by the SRP. The principal 
changes are in the estimator gains. It probably 
would have been possible to obtain a similar sensi- 
tivity reduction by using only the estimator parameters 
as free parameters. The same effect was also observed 
in Example 1. 

C-4 Conclusions 

(a) The method described in this section can provide a consider 
able reduction in the system sensitivity to parameter variations. If 
the initial system is optimal , the output and control rms values will 
increase due to the sensitivity reduction. However, the nature of 
the optimality of the nominal system has to be considered carefully 
since it depends on the assumed values of the process and measure- 
ment noise intensities. In the examples that were examined, the 
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sensitivities were reduced by a factor of 2 to 3, while the output 
rms increased approximately by a factor of 2. 

(b) Both examples given in this section use state estimate 
feedback controllers. In this case the free parameters are the feed- 
back gains and the estimator parameters. The method, however, is by 
no means limited to such controllers. It is equally applicable to 
systems with classical compensation networks or other designs that can 
be represented in state variable form 

(c) The computation time for a 12th order system with 20 
variable parameters was four minutes on an IBM 360/67 computer. 

This computation time is almost insensitive to the number of free 
parameters. It is therefore recommended to define all the 

parameters that are at the designer’s disposition as free parameters, 
at least for preliminary runs. If some parameters do not vary xn 
those runs, they may be fixed for subsequent runs. 

(d) In applying this method, it is recommended to initially 
select the parameter covariance matrix so that the additional Pl(J^) 
is much larger than the nominal Pl(J Q ). A design that is mainly deter- 
mined by sensitivity considerations results. If the nominal properties 
of this design are unsatisfactory, other designs with lower sensitivity 
weightings may be executed. 

The comparison between these designs is made conveniently by 
means of the stability range and the output and control rms values. If 
the system has specific performance requirements such as time response 
envelope or gain and phase margin, these can also be used as comparison 

criteria. 

The most satisfactory design, in general, is a matter of subjective 
designer preference. 


V. DESIGN OF A CONTROLLER FOR A TRACKING TELESCOPE 


A. INTRODUCTION 

In this Chapter, the design of a controller for a tracking tele- 
scope is described. It was selected to represent a commonly encountered 
design problem and is a good example for the use of the design methods 
of Chapter II. The data and the requirements were obtained from various 
sources, mainly by private communication, and typify actual systems al- 
though they do not necessarily represent a specific one. 

Various aspects of the problem of precision pointing and tracking 
have been treated in the literature and several system designs have been 
described [CA-1, FI-1, JO-1, WH-l]. The detailed design of a controller 
for such a system has to consider many specific aspects of the system 
such as structural details, actuator and sensor characteristics, disturb- 
ance inputs, etc. If these aspects are not considered, the system per- 
formance will be degraded, especially for a mobile system which operates 
in a severe and/or changing environment. The design problem is therefore 
involved and is, in great part, of computational and experimental nature. 

Since the principal purpose of this chapter is to demonstrate the 
application of the design methods of Chapter II, a detailed design for a 
changing environment is beyond its scope. Even for a simplified ground- 
based system, structural considerations must be neglected. The actuator 
and sensor characteristics are, however, taken into account and a plausi- 
ble disturbance input is assumed. The design of the controller for this 
simplified but non- trivial system is less involved but the design methods 
can easily be extended to more complicated tracking problems. The sim- 
plified plant is described in Section V-B. The controller specifications 
are defined in Section V-C, and the controller design is described in 
Section V-D. The performance of the controller is examined in Section 
V-E. ' 
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B. PLANT REPRESENTATION 

B-l Description 

The plant is a ground based telescope, the purpose of which is 
to track a moving object (target) in the sky. For angular freedom of 
motion, the telescope is mounted in a three-gimbal structure: inner 

azimuth, elevation, and outer azimuth. It has unlimited angular free- 
dom in azimuth and -6° to +25° freedom in elevation. A schematic 
view of the telescope is shown in Fig. V-l. 



cp inner azimuth gimbal 
angle (limited to 3°) 

S elevation gimbal angle 
(limited to -6° to 25°) 

\|f outer azimuth gimbal 
angle 



FIG. V-l TRACKING TELESCOPE ANGLES 
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Torquers are mounter: on all three gimbals. The inner azimuth 
and the elevation axes have hydraulic torquers. The outer azimuth axis 
has an electric torquer. The transfer functions from torque command 
to torque output for the hydraulic torquers are derived in Appendix C. 
The transfer function of the electric torquer is a pure gain. Since 
the torquers are subject to saturation acceleration, limiting networks 
are mounted at their inputs . 


The purpose of the two azimuth gimbals is to enable the inner 
azimuth torquer to operate on a lower inertia and thereby to provide the 
higher accelerations required for high bandwidth tracking. The outer 
azimuth gimbal then provides mainly the slewing capability. The relative 
angular freedom of the inner azimuth gimbal is + 3°. 


The following measurements are available: (i) the components in 
azimuth and elevation of the error angle between the target line of 
sight and the optical axis of the telescope are measured by detectors 
mounted in parallel to the telescope axis. These measurements are 
sampled at a rate of 120 meas/sec and are held in a zero order hold. 

(ii) The integrals of the angular rates about the inner azimuth and 
elevation axes are measured by rate integrating gyros with their input 
axes along these directions. (iii) The relative gimbal angles are meas- 
ured by resolvers. Numerical data for the system are given in 
Appendix D. 


B-2 State Representation 

The dynamic equations are .derived in Appendix F. Defining 

• ■ • f 

a = cp + \|r cos e (the component of the total angular rate about the lb 
axis) 




+ I^/Oi sin e + ( 1^ - 1^ - sin £ cos e 

2 2 

[l 2 sin e +(l g - I 1 )cos e + + 1^ cos e 

+ ii, + 1., - i,:-h sin e cos € - I <3je cos e 

X. ct O J_ 


(5.1a) 


(5.1b) 

(5.1c) 
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where M , M , and M are the external torques acting about the inner 
a' \|r 

azimuth, elevation, and outer azimuth axes. These torques consist of. 
the torquer outputs and disturbance torques. The moments of inertia 
are defined in Appendix D. 

To linearize Eq. (5.1), the following assumptions are made: 

(a) The elevation angle remains small so that sin £ « e, cos £ 1. 

(b) The terms that consist of products of two angular rates and sin £ 
may be neglected. This is equivalent to assuming small deflections of 
the gimbals from their required orientations. 

The coefficient of \jr in Eq. (5.1c) can be written as 

I 2 sin 2 £ + ( I 3 - I 1 ) cos 2 £ + I 4 

- *3 - h + *4 + (I 2 ' r 3 + I l )si " 2 e 
= I +A1 sin 2 £ - 2750 + 160 sin 2 c 

(from the data of App. E). 


The linearized equations 

are 


M 

= i,5 

(5.2a) 

a 

i 

M 

T *' 

= Le 

(5.2b) 

e 

3 


M 

= V + M a • 

(5.2c) 


To put these equations in state form, the torquer equation for the inner 
azimuth and elevation gimbals as derived in Appendix C, have to be added 
Considering Eqs. (C-8) and (C-9) of Appendix C, it can be seen that the 
linearized elevation equation (5.2b) is decoupled from the azimuth 
equations. The elevation control, therefore, can be treated separately 
from the azimuth control. Using Eqs. (5.2a) and (5.2c) and Eq. (C-8) of 
Appendix C, the azimuth state equations can be written as 
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where (see App. D) 





V 




and w , w and w are disturbance torques and actuator noises. 

In this linearized system, a is the total azimuth angle of the 
inner gimbal, measured from some reference direction. To add the output 
equation, (5.3) has to be augmented. The azimuth error detector measures 
the angle CL - CL , after it is sampled and held in a zero order hold. 
Such a sample-and-hold operation introduces phase lag, and if it is not 
taken into consideration in the design, a system with insufficient phase 
margin may result. A linearization of the sampler and hold is done in 
Appendix G. The linearized sampler and hold is represented by an addi- 
tional state equation 


p = -4f P - 7,2f (a - a ) (5.4) 

o o c 

where f is the sampling rate (120 samples per second). 

o ■ ■ . 

The measurement is given by 




y 


( 5 . 5 ) 


= p + o.9(a - a ) . 

The effect of using the linearized rather than the exact representation 
of the sampler is checked by simulation. Replacing the state q by 


r = a in order to get fewer parameters in the state equation, the final 
state equation (without disturbances) is obtained as 



o 
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Substituting the numerical values of the entries in the F matrix (from 
App. C), the open loop eigenvalues (in rad/sec) of the system are found 
to be 

X = -80 + 615 j (torquer) 

1,2 - 

Xg = -480 (sampler-hold) 

X, = - 0.024 

4 

le a n ~ 0 • 

5,6,7 

The parameter b^ in Eq. 5.6 represents the effect of the outer gimbal 
motion on the inner torquer output. If this term is neglected, X 4 = 0 
results. Relative to the other system eigenvalues, the change in X 4 
is small and the term b therefore can be neglected. This eliminates 

O 

the coupling from the outer azimuth gimbal into the inner azimuth 
gimbal and decouples the azimuth system into two single input systems, 
the controllers for which can be designed separately. 

By linearizing and neglecting small terms, the coupled 
system described by Eq. (5.1) has thus been decomposed into three 
single input subsystems. The inner azimuth and elevation subsystems 
are almost identical and the controller that is designed for one of 
them will also be suitable for the other one with small numerical changes. 


The inner azimuth system is shown below. 



A 
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o 

Oj 

X - 

’.9 

a + 

v d 

_0 

10 0 0 


0_ 

C 

V 
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Structure of process and measurement noise intensity matrices: 



R = 


r d 0 


Numerical values: 


m. 


in = 


160 sec 


-1 


-2 


<i n = 


385,000 sec 

120 samples/sec 

-4 2/ 3 

10 rad /sec 

. -5 '2/ 3 

10 rad /sec 

-11 2 , 

10 rad /sec 


from App. C 


•\ 


\ from App. D 


r = 


“14 2 

2 X 10 rad /sec . J 


Transfer functions: 


0!(s) 

V s> 


a(s) 

u(s) 


2 2 ■ 

s (s + m^s + m 2 > 


(5.8) 


(5.9a) 
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(5.9b) 


g(s) _ i 

W (s) “ 2 

1 s 

The controller for the outer azimuth gimbal has relatively low 
performance requirements and its design is straightforward. 

The linearization and simplification of complex systems as shown 
in this section is common in the design of controllers. It should be 
attempted, whenever possible, since it simplifies the design of the 
controllers to a large extent. It is important to bear in mind, however, 
that as a final design step, a simulation of the full nonlinear system 
should be made using the designed controllers. This simulation should 
be as close as possible to the real system, including all the non- 
linearities and neglected terms. Due to the limited scope of this 
Chapter such a simulation will not be performed here. 

C. CONTROLLER SPECIFICATIONS 

For a meaningful evaluation of the controller design, some quan- 
titative specifications are required. These specifications should 
reflect performance requirements such as tracking and disturbance re- 
jection, as veil as hardware limitations such as component saturation. The 
performance requirements in this Section are based on typical requirements 
of actual tracking systems. For such systems; tracking requirements are 
often defined in terms of the permissible error when tracking a target 
moving in a straight line [FI-1 j. It is also customary to require that 
no steady state output error occur when constant disturbances such as 
steady wind or unbalance torques act on the system. 

The numerical values given in the specifications reflect levels 
of performance that are plausible with the assumed process and measure- 
ment noises, and with the given detector sampling rate. The influence 
of different noise level assumptions are discussed later. According 
to this approach, the controller specifications are formulated as follows, 
(a) Tracking. The system can track a target moving with constant angular 
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acceleration with finite steady state error. It is desirable for this 
error to be less than 100 [trad, which is the value of the detector 
rms noise, for the entire acceleration range of the outer gimbal 
(0.65 rad/sec 2 ). This error corresponds to an acceleration error co- 
efficient of 150 (lrad/( rad/sec 2 ) . This requirement is a simplification 
of the tracking requirement for a target on a straight line flyby. 

The connection between these two requirements is derived in Appendix 

E. 

(b) Torquer saturation . The torquer input is less than 50% of 
the acceleration limit when no tracking is required (pointing at a 
stationary target). This specification is required in order to avoid 
saturation of the torquer by disturbances and measurement noise and 
to retain sufficient control authority for the tracking. 

(c) Large signal operation. The system remains stable for large 
output errors up to 30 mrad. This requirement arises because the tor- 
quers are subject to saturation (V-B-l). At large error signals, 

such as may arise during acquisition, saturation instability therefore 
may occur (Section V-E). 

(d) Steady disturbance rejection . No steady state pointing error 
occurs when steady disturbance torques are acting on the system. The 
maximum transient error for such a disturbance is 

— < 100 iirad/( rad/sec 2 ), 

8w 

where w may be either or w g (Eq. 5.7). With this numerical value 

a disturbance torque of 25% of the inner azimuth torquer capacity 

will cause a maximum deflection that is of the order of the measurement 

noise. 

’ (a) Transient response . The system has "well behaved" tracking 

and disturbance transient responses by common engineering criteria. 

The tracking bandwidth is limited by the sampling rate. A bandwidth 
of 25% to 50% of this rate is considered reasonable. 
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D. DESIGN OF A CONTROLLER FOR THE INNER AZIMUTH 


GIMBAL 


Dt- 1 Introduction 

The design procedure for a state estimate feedback (SEF) con- 
troller (Ch. II) consists of the following steps: (a) estimator 

design, (b) controller design, (c) evaluation* In general, several 
systems are designed in the first two steps, each one having different 
state and control weights and, in some cases, different estimator gains. 
The nature of the evaluation depends on the system specifications. 

In general, it consists of several stages, and in each stage, some de- 
signs are eliminated. In some cases, more than one design cycle is 
required in order to satisfy the system requirements. 

The three-step design procedure outlined above is fairly standard. 
The nature of each one of these steps, however, depends on the system 
specifications and on the special features of each specific systems. 

In the inner gimbal system, the following special features have to be 
considered: (a) although the system is SISO it has two measurements 

and different estimators can therefore be designed using various combi- 
nations of those measurements. The various designs are described in 
Section V-D-2. (b) The system has to be augmented by integral control 

in order to satisfy the tracking and disturbance rejection require- 
ments (requirements a and d in Sec. V-3). The need for this aug- 
mentation is established in Section V-D-2 and the controller for the 
augmented system is designed in Section V-D-3. (c) The evaluation has 

to include the determination of stability for large input . This is 
done in Section E-l .. . : 

In Section D-4, the various possible designs are compared accord- 
ing to the linear performance criteria and two designs are selected 
for sensitivity comparison. This comparison is made in Section D-5 
and as a result, a final system for further performance evaluation is 
selected. This evaluation is described in Section E-5. 
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D-2 Controller and Estimator Structure 


Two measurements are available in the system: (a) target detector, 

(b) gyroscope. The gyroscope can be used in one of two modes: as 

a rate integrating gyro, or as a rate gyro. The dynamics of the gyro- 
scope in its two modes are described in Appendix H, 

The performance specifications call for (a) finite steady state 
output error for constant acceleration command; (b) zero steady state 
output for constant disturbances. The satisfaction of these require- 
ments depends on the structure of the system and not on the values of 
its parameters. The level of the steady state acceleration error and 
the maximum transient deflection for constant disturbance depend on the 
parameters . 

Several estimator designs will be discussed in this Section, and 
the structures of the systems incorporating them are examined. The 
numerical performance criteria such as output and control noise, 
maximum disturbance deflection, and steady state acceleration error 
can only be evaluated after both the controller and the estimator are 
designed. This evaluation is therefore deferred to Section D-4. The 
discussion in this Section is limited to estimators, the order of which 
is at most, equal to the order of the system. 

( 1 ) Estimator using the detector measurement only . Equations 
(2.58) and ( 2 . 67 ) define the order of the polynomial equilibrium control 
required for the tracking of polynomial reference outputs and rejection 
of polynomial disturbances. Using these equations on the transfer 
functions of Eq. (5.9) and referring to the controller specifications 
(V-c) , the results are: (a) a constant acceleration reference output, 

can be tracked with constant equilibrium control. No integration of 
the output error is therefore required in order to obtain a constant 
output error for this reference output. (b) A constant disturbance 
requires a constant equilibrium control for zero output error. One 
integration of this error is therefore required. The disturbance re- 
jection specification, therefore, determines the requirement for in- 
tegral control. With this control the output error is zero for a 
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constant acceleration reference output. In addition to the equili- 
brium control, error feedback control is required. From Chapter II-D-3 
the system with this reference output is completely reducible, i.e. , 
the error state has the same dynamics as the system state. It can 
therefore be estimated by an estimator that uses the model of the 
system, the output of which is compared to the target detector output. 
This estimator can be used for generating the error feedback gains. 

For zero target motion, it obviously generates the state estimates. The 
state equations of this estimator are given below. 



The structure of a system using it is shown in Fig. V-2. 


A controller using this estimator more than satisfies the steady 
state error requirements since it has zero acceleration error where only 
a finite error is required. Its transient disturbance response may, 
however, not be satisfactory. The detector noise level is relatively 
high (see App. D) and the estimator bandwidth is therefore low. 
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FIG. V-2 SYSTEM WITH DETECTOR MEASUREMENT ONLY 


The closed loop transfer functions from the disturbances to the 
output contain in their denominators, the eigenvalues of both the con- 
troller and the estimate error systems. The transient response to these 
disturbances may therefore be dominated by the estimate error eigen- 
values and therefore may be large and relatively slow. In order to 
improve this response, the gyroscope can be used as an additional meas- 
urement in one of its two modes. Estimators using this additional meas- 
urement are described below. 

(2) Estimator using rate integrating gyro as measurement . As ex- 
plained in (1), a system using an estimator with detector measurement 
only may not have satisfactory disturbance response. This is so be- 
cause the information about output errors due to disturbances is obtained 
through the detector measurement only, which is relatively noisy and 
therefore must be filtered by a low bandwidth filter. The gyro has a 
much lower measurement noise than the detector (App. D). 

The disturbance information obtained when the gyro is used as 
an additional measurement may therefore be filtered through a faster 
filter and higher stiffness to disturbance outputs may thus be obtained. 

The output of a rate integrating gyro is an integral of the torques 
acting about its Output axis. They are: (a) torquer generated torques, 

T, ( b) drift torques, D, (c) torques caused by angular rates about the 
input axis, coH. Therefore, 
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cp 


k D (T + D + uH) . 


(5 . 10 > 


The total output must be kept close to zero. A command must therefore 
be applied to the torquer such that it generates torques that balance 
the drift torques and those due to angular rates. This command may be 
obtained from the detector signal e such that 

f (e) + D + oH = 0 (5.11) 

where f(e) is an unspecified function of the signal e. Since, for 
constant angular rate, a zero error is required, this function must 
include the integral of the error and can be of the form: 

f(e) - a Q e + | e dt . (5.12) 

The gyro-detector combination is shown in Fig. V-3. If the 
gyro is used as the sole measurement for the estimator, the gyro output 
will be kept close to zero by the controller, which adjusts the value 
of e so that the gyro torques are balanced. The single gyro measure- 
ment can be considered to measure a linear combination of system states 
only if the target motion is viewed as noise, which is a reasonable 
assumption for an unknown maneuvering target. In that case the measure- 
ment is a linear combination of the state a and its two integrals. In 
order to model this measurement correctly in the estimator, the system 
model has to be augmented by the two integrals. A seventh-order estim- 
ator thus results. Since the order of the estimator is required not 
to exceed the order of the system, such an estimator cannot be used. 

An estimator of the order of the system is obtained by using the model 
of £q. (5.7) and considering the gyro as a measure of the state a. 

Since in reality this measurement also contains the error integrals, the 
estimator does not estimate the actual states and a coupled system re- 
sults. The block diagram of this coupled system is shown in Fig. V-4. 
From this figure, it can be seen that the effect of the non-modelled 
integration is to close an outer loop around the controller from the 
zero order hold output to the gyro torquer input. Its transfer function 
is: 
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This transfer function includes one integration in the gyro. Note that 
the open loop transfer function of this augmented system has two roots 
at the origin only after integral control has been added. Therefore, 
integral control is required in order to obtain finite steady state 
acceleration error. This is to be compared with the system in which the 
estimator uses the detector measurement only. There, zero steady state 
acceleration error is obtained with integral control. 

It is to be noted that because of the additional measurement, the 
steady state tracking and disturbance rejection properties are not 
obtained from the open loop transfer function of the actuator and plant. 
From Fig. V-4 it can be seen that the additional measurement has the 
effect of closing an inner feedback loop. The closed loop transfer 
function of this inner system is considered as the open loop t. f. for 
the determination of the steady state properties. If the constants 
a Q and a_ L are small (a Q /a 1 much smaller than the eigenvalues of the 
nominal controller), the outer loop may be neglected and the feedback 
and estimator gains found by OPTSYS. The actual system eigenvalues 
will be somewhat shifted from their assumed values. If, however, rela- 
tively tight integral control is required in order to obtain a suffi- 
ciently small acceleration error, the feedback and estimator gains and 
the constants and a^ have to be found by parameter optimization. 

The augmented system shown in Fig,, V-5 is used for this optimization. 

Is is performed using the program PAROPT (Ch. IV-B-4) without the 
sensitivity reduction option. 

The weighting matrices and the covariance matrices are the same 
ones that are used in OPTSYS, 

Although the parameter optimization method is somewhat less con- 
venient to use than the optimal controller and filter method, chis, 
by itself, is not a serious drawback of this sytem since the optimization 
is only carried out a relatively small number of times. 
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FIG. V— 5 SYSTEM WITH ESTIMATOR USING RATE INTEGRATING GYRO OUTPUT EFFECT 

OF NONMODELED INTEGRATION. 


(3) Estimator with rate gyro as sensor . The transformation of a 
rate integrating gyro into a rate gyro by caging it, i.e., feeding back 
its output to its torquer, is described in Appendix H. The rate gyro 
is a second order system for which the natural frequency is determined 
by the caging loop gain and the damping by the rate integrating gyro 
time constant. If the natural frequency is sufficiently higher than 
the system natural frequencies, the gyro transfer function can be 
approximated by a constant, viz., 


cp « k w . (5.13) 

r cr 

The effect of neglecting the gyro dynamics in the estimator model is 
to couple the estimator and system eigenvalues. The amount of coupling 
depends on the natural frequency of the gyro. To minimize the coupling, 
this frequency has to be high but this increases the gyro measurement 
noise. In Appendix H these effects are discussed and a natural frequency 
is selected. The conflict between high noise and coupling can be avoided 
if the rate gyro is modeled in the estimator but this leads to a more 
complicated estimator (7 states). 

Using the approximation of Eq. (5.13), an estimator can be de- 
signed using both the detector and rate gyro as measurements. As ex- 
plained above, the addition of the rate gyro should improve the disturb- 
ance response. 

Three different realizations of this estimator are considered. 

(a) One full state estimator using the two measurements. The 
block diagram of this estimator is shown in Fig. V-6a. (b) One reduced 

order estimator using the same measurements. Its block diagram is shown 
in Fig. V-6b. Since its design method is substantially different from 
the one used for the other estimators, it is described in Appendix K. 

(c) Two separate full state estimators: • one two-state estimator for 

q: and (3 (see Fi. V-3a) using the detector as the measurement and the 
rate gyro a.s a known input; • one three-state estimator for w, a, and 
r using the rate gyro as a measurement. 
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FIG. V-6a FULL STATE ESTIMATOR WITH TWO MEASUREMENTS 



FIG. V-6b REDUCED ORDER ESTIMATOR WITH TWO MEASUREMENTS 




















The block diagram of this estimator is shown in Fig. V-6c, 
dynamic equations of its models are given in Fig. V-7. 
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FIG. V. 7a MODEL FOR ESTIMATOR NO. 1 
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FIG. V-7b MODEL FOR ESTIMATOR NO. 2 
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The 


Because the system is SISO ; the steady state tracking and disturbance 
rejection properties can be studied from the equivalent transfer func- 
tion representation of the estimators. 

All three designs can be represented by the same transfer function. 
The block diagrams of this transfer function are shown in Figs. V-8a and 
V-8b, where Fig. V-8b is a simplification of Fig. V-8a. For the full 
state estimators, = 0 and = 0. 



FIG. V-8b 


FIGS. V-8a, b EQUIVALENT BLOCK DIAGRAMS OF SYSTEM WITH RATE 
GYRO AND DETECTOR MEASUREMENT. 
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From Fig. V-8b, it is observed that the equivalent open loop 
transfer function, which contains the effect of the integral control 
and the rate measurement feedback, has two roots at the origin. The 
systems can therefore track a constant acceleration command with constant 
control and finite tracking error. To obtain zero tracking error for 
acceleration command, one additional integration has to be added. 

The effect of the gyro measurement, both when used as a rate integrating 
gyro or a rate gyro, is therefore to reduce the system from type (2) 
to type (1). This disadvantage will be shown to be outweighted by the 
improved transient response and disturbance rejection discussed in 
Section D-4 (below). 

The five designs that were described in this section cover all 
the possible combinations of the measurements. There may, however, 
exist other designs using these same measurements. The structure of 
the described designsis shown in Table V-l, Later, the gains are cal- 
culated and the designs compared. 

Table V-l 


ESTIMATOR STRUCTURES 


DESIGNATION 

MEASUREMENTS 

- . . 

METHOD OF DESIGN 

j|g[|j| 

REMARKS 

D 

detector 

Separate design of 
controller and | 

estimator by OPTSi’5 

1 

1 

2 ! o 


D1 

detector 
and rate 
integrating 
gyro 

Combined design by 
parameter optimi- 
zation 

m 

Coupling between 
control ler and 
estimator roots 

m l 

detector and 
rate gyro 

same as type D 

1 j 0 


DR 2 

detector and 
rate gyro 

same as type D 

i i o 

two separate 
estimators 

DRR 

detector and 
rate gyro 

parameter optimi- 
zation (see App. J) 

i I 0 
1 

reduced order 
estimator 
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D-3 State and Control Weight Selection 


In this Section the determination of the state and control weights 
is described. These weights are then used in an optimal control program 
for the determination of the feedback gains and closed loop eigenvalues 
and eigenvectors. Several systems are designed for subsequent evalua- 
tion in the next section. 

The relation between the weights and the system eigenvalues is 
determined by the root square locus method. This method was mentioned 
in Chapter II and is described in more detail in Appendix It en- 

ables one to determine the weights fox' single input systems so as to 
obtain required eigenvalues. The required eigenvalues, however, have 
to be determined from the given system specifications. The corres- 
pondence between these specifications and the eigenvalues is not unique 
since the time response also depends on the controller eigenvectors and 
on the estimate error eigensystem. A rough correspondence can, however, 
be determined by assuming that the controller and not the estimator 
dominates the time response. Under this assumption, the following 
criteria for the eigenvalue locations can be established. (n) n high 
bandwidth is desirable for low tracking error and good disturbance re- 
jection. It is limited by the detector sampling rate (W g = 800 rad/sec) 
and the actuator noise requirements. A bandwidth of w = 100 sec 
seems reasonable as an initial point. (b) The actuator has a natural 
damping coefficient of 0.13. Although this is somewhat low, high 
gains may be required in order to increase it and therefore initially 
this location will be considered satisfactory. 

These two criteria sum up to the requirement that the damping 
of the actuator 1 roots should not decrease and that the magnitude of the 
real part of the additional roots should be greater than 100 sec . 

In the remainder of this section, the selection of weights for several 
systems will be described and candidate systems for the performance 
evaluation will be determined. 

The open loop system for which the controller is designed is the 
system of Eq. (5.7) augmented by an integral state i such that 
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- y = p + 0.9a . 

at 


The requirement for this state was explained in Section D-2. In order 
to determine which states have to be weighted, the following rule is 
used [ AN- 2] ; to ensure stability of the closed loop system, the pair [F,E>] 
has to be observable, where D is any matrix such that DD = A, 

From this rule it is obvious that the only state that must be weighted 
is the integral state i. The effects of the weighting of the states 
01 and to are also examined. A multiple-parameter root square locus is 
constructed by first constructing the root square locus for the state 
i , selecting a temporary weight for this state and considering the closed 
loop eigenvalues obtained for this weight as denominator eigenvalues 
for the a root square locus. The same procedure is followed for obtain- 
ing the to root square locus. 


The transfer functions to the three weighted states are: 


G_^ (s ) 


i (s ) 
u(s) 


0. 9m^ (s - 4f Q ) 


N. (s) 

A 1 

~3 ~2 . " D (s) 

s (s + 4f )(s 4- mjS + m 2 ) 


(5.14a) 


V £) 


, x m s(s + 4f.) 
q(s) _ 2 o_ 

u(s) " D.(s) 

i 


(5.14b) 


G to ( s) 


ms (s + 4f ) 
co(s) _ 2 0_ 

u(s) “ D a (s) 


(5.14c) 


where D. (s) is the characteristic polynomial obtained when the integral 
state only is weighted with the selected weight, and s) is the 

characteristic polynomial obtained by weighting both i and d. The root 
square locus for the i, d, and to states are shown in Figs. V-9 through 
V-ll. From these figures, the influence of increasing the weights 
of the different states can be assessed as follows: 
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FIG. V-9 ROOT SQUARE LOCUS FOR INTEGRAL STATE WEIGHTING. 


Integral Weight: = 2X10 sec 



FIG. V-10 ROOT SQUARE LOCUS FOR POSITION STATE 
WITH FIXED INTEGRAL WEIGHT. 
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FIG. V-Xl ROOT SQUARE LOCUS FOR RATE STATE WITH FIXED INTEGRAL AND 
POSITION WEIGHTS. 



(a) Integral state: • little effect on actuator roots; • other 

roots move away from the origin with little change in damping of com- 
plex roots, (b) Position: • some effect on actuator roots; • real 

root approaches origin; • complex roots move away from the origin; 
damping of the complex roots increases; (c) Rate: • increased 

bandwidth and damping of actuator roots; • real root moves away from 
the origin; • complex roots initially remain at same distance from 
imaginary axis with increased damping. 

Roughly, therefore, the integral state weight influences primarily 
the bandwidth, whereas the rate weight affects mainly the damping. The 
position weight has an intermediate effect. Using these considerations, 
three designs were selected for further evaluation as shown in Table 
V-2 . 


Table V-2 

CONTROLLER DESIGNS 


Design 

No. 

Characteristics 

(sec -6 ) 

0 

a 

(sec -4 ) 

a w 

(sec - 2) 

Remarks 

i 

nominal (N) 

13 

2 X 10 

0 

0 

Point 2 in Fig. V-9 

2 

high gqin (HG) 

2 X 10 14 

0 

0 

Point 3 in Fig. V-9 

3 

high damping 
(HD) 

2 X 10 13 

2 X 10 9 

4 

8 X 10 

Point 2 in Fig. V-ll 


The selection of Design 1 as an initial design was made by considering 
its acceleration error. If state feedback, without an estimator is used 
the steady state acceleration error is: 

e rate gain 

a ~ integral state gain 

c 

For Design 1, this ratio is 1.2 x 10~ 4 sec 2 , or 120 (j,rad/( rad/sec ) . 
This is less than the required acceleration error but an increase in 
■this Grror is to b6 0 xp Get © d if an. cstitn&tcir is used "fco obtain the 
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feedback gains (sea App. K). 

Design 2 has a lower acceleration error, Design 3 a higher one. 

No design using integral and position weights only was evaluated since 
this represents an intermediate case between the cases that were 
selected . 

D-4 Estimator Design Evaluatio n. 

The following designs were made and are evaluated in here. 

(a) The controllers and estimators given in Tables V-l and V-2 
with the noises as given in App, D and H. The influence of different 
noise assumption is discussed in Section E-3. (b) An additional 

estimator of type DR 2 (two estimators, one using the rate gyro meas- 
urement, and the other the detector measurement). In this design 
(DR 2H) , the gains of the detector estimator were arbitrarily selected 
at higher values in order to decrease the acceleration, error. The 
connection between the estimator gains and the acceleration error is 
developed in Appendix K. 

The gains for the controllers and the estimators (except Type 
DI and DRR) were found by the optimal control program, OPTSYS [ BRY-3] . 
The feedback and estimator gains for Type DI and the estimator gains 
for Type DRR were found by parameter optimization using the sensitivity 
minimization program, PAROPT, described in Ch. IV-B-4, without the 
sensitivity reduction option, and with the same state weights as were 
used for Design 1 in Table V-2. The eigenvalues of the controllers 
are shown in Figs, V-9 to V-ll. The eigenvalues of the estimators 
are shown in Fig. V-12a, The eigenvalues of the coupled system DI are 
shown in Fig, V-12b. It is observed in Fig. V-12a that the reduced 
order estimator has a very slow root. This can be made plausible by 
considering Eq. J-15 in App. J . Comparing this equation to Eq. (2.31), 
it can be seen that in the full state estimator, the measurement noise 
forces the estimate error only, whereas in the reduced order estimator, 
this xioise forces both the state and the estimate error through the 
estimator gain K. For low output and control noises, a lower gain 
K is therefore desirable and slower estimator roots result. 
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FIG. V-12a ESTIMATE ERROR SYSTEM EIGENVALUES FOR DIFFERENT 
TRACKER DESIGNS ■ 
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The performance of various combinations of these designs was 
evaluated according to three criteria: (a) steady state error for 

acceleration input; (b) maximum deflection for step disturbance 
input; (c) control mean square noise. 

For systems using the DI and DR 2 estimators, it is relatively easy 
to find the analytic expression for the acceleration error. This is done 
in Appendix K. For systems using the DR 1 estimator, the analytic 
expression is complicated and the error was calculated numerically 
using the design program XAGSA [WIT-1] * For systems using the Type D 
estimator, there is no steady state acceleration error. 

The disturbance deflection was found in all cases with the aid 
of the XAGSA program. The control rms values were found with the aid 
of the OPTSYS program. The results of these evaluations are given in 
Table V-3. The time responses of some of the systems are shown in Figs. 
V-13 and V-14. 


Table V-3 

PERFORMANCE COMPARISON OF CONTROLLER AND ESTIMATOR DESIGNS 


Design 

No. 

Estimator 

Type 

Controller 

Type 


Maximum De- 
flection for 
Step Disturb. 

( 4 rad v 

Control 
rms noise 

(rad/sec^ ) 

Output 
rms noise 

(grad ) 


\ rad/ sec^ ) 

a 

D 

N 

0 

260 

1.8 

65 

2 

DR 1 

N 

280 

145 

0.54 

111 - .!. \ 

29 

3 

■ 

DR 1 

: ■ 

KG 

205 

96 

.V • ■' ■ ' ; 

0.94 

26 

4 

DR 1 

HD 

590 

180 

0.58 

30 

5 .. 

DR 2 

N 

400 

48 - 

1.05 


6 

DR 2H 

N 

185 

48 

1.8 

46 

m 

' DR 2H 

HG 

130 

not Computed 

3.3 

53 

8 

DRR 

N 

28000J : 

18 

BSH 

28 

9 

DI 

N 

240 

not computed 

2.5 

78 


N = nominal HG = high gain HD = high damping 
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FIG. V-13 OUTPUT ERROR FOR CONSTANT ACCELERATION COMMAND 




e 



FIG. V-14 CONSTANT DISTURBANCE RESPONSE OF DIFFERENT 
TRACKER DESIGNS. 
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A control rms higher than 2 rad/sec 2 i s considered unacceptable 
since it violates the torquer saturation requirement. Other criteria 
for such designs were therefore not evaluated. 

The following observations can be made from Table V-3 and Figs. 

V-13 and V-14. (1) the addition of the rate measurement to the de- 
tector measurement lowers both the control rms and the disturbance de- 
flection (compare Designs 1 and 2), However, it causes the system to 
have a finite, instead of a zero, acceleration error. (2) Increasing 
the bandwidth of the system lowers acceleration error and disturbance 
deflection but increases control rms (compared Designs 2 and 3). This 
result is to be expected. (3) Increasing the damping increases the 
acceleration error and the disturbance deflection (compared Designs 
3 and 4). The time response of the better damped system has no over- 
shoot (see Fig. V-13). (4) Two partial estimators instead of one full 

estimator (with the same controller) cause the system to have higher 
acceleration error and control rms but lower disturbance deflection 
(compare Designs 2 and 5). (5) The designs using the rate gyro as a 

second measurement have generally better performance than the one using 
the rate integrating gyro (DI). (6) The reduced state observer (DRR) 

has a higher control noise than the full state observer with the same 
controller, as is to be expected. However, it also has an unacceptably 
high acceleration error. Its acceleration response is dominated by the 
low eigenvalue = -0.56 rad/sec. 

The ranking of the designs according to the criteria of Table V-3 
is given in Table V-4. The designs below the dashed lines in this Table 
are considered unacceptable for the respective requirements. Note that 
designs with acceleration errors of up to 200 |lrad/(rad/sec) are con- 

: ' . ■ ■ ' - ■ Q 

sidered acceptable, although in Section V-C an error of 150 (Lrad/(rad/sec~) 
is specified. This is so because this error is considered as a design 
goal only and not as a rigid specification. From Table V-3 it can be 
seen that satisfying this requirement would cause an unacceptably large 
controller noise. 
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Table V-4 


RANKING OF THE DESIGNS 


Acceleration 
Steady State 
Error 

Disturbance 

Deflection 

Control 

rms 

1 

8 

2 

7 

5, 7 

4 

6 

3 

3 



5 

O 



9 

2 

8 

2 

4 

6, 1 

v 

1 


9 

% 

9 

8 


7 


Note that in. the above Table, only Designs 3 and 6 are acceptable 
according to all the criteria. Designs 3 and 6 are comparable and the 
selection between them will be made according to parameter sensitivity 
criteria in the next Section. 

Design 8, which uses the reduced order observer, is clearly un- 
acceptable as is. It, however, merits some more consideration since it 
is simpler than all the other designs. It's acceleration error may be 
decreased by shifting the observer poles but that, of course, will 
increase its control and state noise. It is reasonable to assume that 
for parameter sensitivity reduction, the slow eigenvalue has to be 
moved further from the imaginary axis and therefore the applications 
of the sensitivity reduction program to this system may also improve 
its acceleration response. This is done in the next Section also. 
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D-5 


Sensitivity Reduction and Final Selection 


The sensitivity of Designs 3 and 6 to the variation of the hydraulic 
actuator spring constant and damping coefficient is shown in Table V-5. 

The region of stability is considered as the sensitivity measure. From 
this Table it can be seen that both designs have low sensitivities to 
damping variations. The sensitivity to stiffness reduction of Design 
3 is unacceptable, whereas for Design 6, it is tolerable. 


Table V-5 

STABILITY REGION FOR HYDRAULIC ACTUATOR PARAMETER 

VARIATIONS 




Design 3 
(DR l-HG) 

Design 6 
(DR 2H-N) 

Spring 

Constant 

+ Ak/ k 

Large 

Large 

Range 

- Ak/k 

o 

o 

i 

-0.18 

Damping 

- Ab/b 

Large 

Large 

Coefficient 

Range 

- Ab/b 

-0.7 

-1.1 


The sensitivity reduction method may be applied to Design 3, Since 
this is an optimal system, its state and control noises will increase 
due to this reduction but this is acceptable since the control noise is 
well below the specified level. Design 6 is not optimal since the gains 
of estimator No. 1 were selected arbitrarily. Sensitivity reduction 
may, therefore, increase or decrease its control noise. However, a 
system with lower control noise will most probably have a higher accel- 
eration error. Since both this error and the control noise are close 
to their specified limits, the application of the sensitivity reduction 
method to this system is not feasible. The result of the sensitivity 
reduction for Design 3 are given in Tables V-6 and V-7. 
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Table V-6 


STABILITY RANGE OF DESIGN 3 




Original 

Desensitized 

Spring 

Constant 

+ Ak/k 

Large 

+ 0.31 

Range 

- Ak/k 

t 

o 

o 

>d- 

- 0.22 

Damping 

+ Ab/b 

Large 

Large 

Coefficient 

Range 

- Ab/b 

-0.7 

-0.65 


Table V-7 

NOMINAL PERFORMANCE CRITERIA OF DESIGN 3 


System 

■ 

Steady State 
Error For 
Acc. Input 

Maximum De- 
flection for 
Step Disturb. 

Control rms 
Noise 

Output rms 
Noise 


f prad \ 

— ■ 

/ 2 

(rad/sec ) 

(prad ) 


\ rad/sec^/ 


Original 

205 

96 

0.94 

26 

Desensitized 

250 

28 

2.8 

45 


Comparing these Tables with Tables V-3 and V-5, it is obvious that Design 
6 is preferable by both nominal and sensitivity criteria. 

As explained in the previous Section, the sensitivity reduction method 
is also applied to the reduced order estimator design, No. 8. In this case, 
this method is used essentially as a pole placement method. The results of 
its application are given in Table V-8. From this Table it can be seen that 
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the sensitivity reduction method has, as expected, considerabl y decreased 
the acceleration error for this case. The stability ranges of Design 
6 and of the desensitized reduced order estimator, Design 8A, are compared 
in Table V-9. These two designs can be seen to have comparable prop- 
erties. The nominal properties of Design 6 are slightly better but 
the sensitivity of Design 8A to reduction of the spring constant is 
considerably lower. The higher sensitivity of this design to reduc- 
tion of the damping coefficient is of no importance since the stability 
range is more than adequate. Since this design is also simpler, it is 
selected as the final design. 


Table V-8 

SENSITIVITY REDUCTION OF REDUCED ORDER ESTIMATOR 


System 

Steady State 
Error For 
Acc. Input 

Maximum De- 
flection for 
Step Disturb. 

Control 

rms 

Noise 

j Output 
J rms 
< Noise 

Estimator 

Eigenvalues 

/ grad \ 

/ grad \ 

(rad/sec 

1 

» 

2) (|trad*) 

i — 

( rad/ sec) 


\rad/sec “ £ J 

\ rad/ sec */ 

Original 

28000 

18 

1.2 

t- 1 

| 28 

- 0.54 
-36 4 610j 

Desensitized 

195 

40 

2.1 

| 

! 43 
i 

i 

- 110 

- 94 4 60Gj 


It is important to note that since the system is noise limited, a 
different system may have been selected if different values were assumed 
for the measurement noises and disturbances. In a more detailed design 
procedure, several designs are usually made using different noise levels 
and ohe of the considerations in the selection of the design may also 
be low sensitivity to assumed noises. 

■In Fig. V-15 , the frequency response to a tracking command is 
shown. It has a bandwidth of 190 sec 1 which is about 25% of the 
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Table V-9 


STABILITY RANGE COMPARISON OF TWO ESTIMATOR DESIGNS 




Design 6 

Desensitized 
Reduced Order 
Estimator 

Spring 

Constant 

+ Ak/k 

Large 

Large 

Range 

- Ak/ k 

-0.18 

-0.33 

Damping 

Coefficient 

- Ab/b 

Large 

Large 

Range 

- Ab/b 

-1.0 

-0.7 


sampling frequency. This is within the range of the specifications 
discussed earlier (V-C), The output response to a constant acceleration 
command and to a constant disturbance is shown in Fig. V-16. 

D-6 Summary 

Among the Designs examined in this Section, two have the best 
overall performance, with the assumed measurement and process noises: 

(a) Design 6, DR 2H. A system using two separate estimators 
with the gains for the first estimator determined by pole 
placement . 

(b) Design 8, DRR, desensitized. A system using a reduced order 
estimator. The sensitivity reduction for this design was 
only used as a pole placement device. 

Design 8 was finally selected. Its eigenvalues, performance criteria, 
and sensitivities were given, in Table V-8; its frequency and time 
responses were shown in Figs. V-15 and V-16. 
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E. PERFORMANCE OF THE SELECTED CONTROLLER 


E-l General 

Some aspects of the performance of the selected controller will 
be examined in this Section. Its large signal operation is examined 
in the next Section and a nonlinear compensation network is introduced. 

In Section E-3 the effect of the sampler and hold linearization is 
assessed. In Section E-4, the tracking performance is compared with 
that of a similar system using aided tracking, and the effects of assum- 
ing a lower measurement noise are evaluated. 

E-2 Large Signal Operation 

As described in Section V-B, an acceleration command limiter is 
placed at the input to the torquer. The purpose of this limiter is to 
protect the torquer from saturation. Its characteristics are shown 
in Fig. V-17, It is not taken into account in the linear design of 
Section V-E, the results of which are therefore only valid for condi- 
tions in which the control amplitude is below the limit, i.e., for small 
displacements from the equilibrium condition. During target acquisition, 
however, the error amplitudes may be large and the acceleration command 
may exceed the limit. The stability of the system for large error sig- 
nals, therefore, has to be examined. 
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The block diagram of the system with the acceleration limiter 
is shown in Fig. V-18. Since the control to both the plant and the 
estimator passes through the limiter, the estimator estimates the states 
correctly even when the limiter is operative. This is therefore the 
case described in Chaptn ' II-C-4 of a state estimate feedback single- 
input single-output controller which has a variable control gain. 

Its dynamic equation is 

x = Fx + c^Gu t 

where c^ is a variable scalar. In Appendix B it is shown that such 
a system has an infinite gain margin for c Q , but may become unstable 
if c is decreased. It is also shown that the use of state estimates 
instead of the actual states for feedback has no influence on the sta- 
bility considerations, and the stability analysis can therefore be made 
assuming that state feedback is used. 

For the inner azimuth system, it can be seen from Fig. V-17 that 
for low control amplitudes c Q = 1 and that it decreases with increasing 
amplitudes. Instability is therefore possible. Assuming that state 
feedback is used, the closed loop characteristic equation can be written 
in the form 

3 2 

s (s + s Q )(s + m^s + m 2 > 

4 3 2 (5 - 15) 

4- cm [(s+4f ) (s c + s c +s o + sc + c.) 

02 o p a Co 0! l 

- 7.2f (sc + c . ) ] = 0 

o y l 

from which a root locus as a function of c^ can be constructed . 

From this equation .is is obvious that for small values of c Q , 
the system will become unstable since the open loop system (c^ = 0) 
has three roots at the origin. This is true for any set of the feed- 
back gains as long as integral control is used. Without integral con- 
trol, the open loop system has only two roots at the origin and there- 
fore will be stable for all values of c Q . 


i 
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The Cq root locus for two values of the integral control gain 
is shown in Fig. V-19. The reduction of the integral gain can be seen 
to increase the stability region considerably. Therefore, although the 
large signal instability cannot be eliminated if integral control is used 
it can be made to occur at higher signal levels if the integral control 
gain is reduced. It is therefore desirable to make the integral gain a 
function of the control amplitude so that for low signals, the tight 
integral control required for low tracking error is preserved. This is 
hard to implement in an analog system but if the control signal is 
assumed to be roughly proportional to the error signal, the same effect 
can be obtained by placing a limiter at the input to the integral con- 
trol integrator. This is shown in Fig. V-20. The characteristics of 
the limiter are: 


y out 

= y in 

for 

y in 

< 

y i 

y out 

II 

for 

y in 


y £ * 


The value of y^ has to be set higher than the error level expected 
for the maximum acceleration command so as not to interfere with the 
acceleration tracking capabilities. The error signal level up to which 
the system remains stable has to be determined experimentally. 

E- 3 Effect of Sampler Linearization 

The state representation of the system as derived in Section B-2 
contains a linearization of the sampler and hold. The effect of this 
linearization has to be determined* The block diagram of V-6c is shown 
in simplified form in Fig. V-21. The dynamics of the continuous system 
can be represented as 

* = V + G c y ’ x <°> = *0 (5.16a) 

where y is considered as a control. The discrete equivalent of Eq, 

<5. 16a) is 
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(a) NOMINAL INTEGRAL GAIN 


jw (sec ) 




FIG. V-19 ROOT LOCUS AS A FUNCTION OF c 


LIMITER 



FIG. V-20 LIMITER FOR THE INTEGRAL GAIN 



FIG. V-21 SIMPLIFIED REPRESENTATION OF CONTINUOUS SYSTEM 
WITH SAMPLER AND ZERO ORDER HOLD (ZOH). 
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(5.16b) 


x. = <t>x. + Ty. , x(0) = x 

i+l ix 0 


where <J> is the transition matrix of the continuous system. Using 


y. = a . - a. = a . - Hx. 

1 Cl 1 CX X 


(5.17) 


Eq. (5.16b) becomes 

x . = ( $ - PH)x. + TO ., x(0) = x n . (5.18) 

x+1 x ex 0 

From this equation the time response can be found for d . and x^. 

ci 0 

This was done with the aid of the computer program SIMUL [ KA-l] 
that determines the state transition matrix $ and discrete control 
distribution matrix, P, from the dynamics matrices F and G. The 
same program also calculates the time response of the system. 

In Fig. V-22, the true time response, computed by this method, 
is compared with the one obtained from the system with the linearized 
sampler and hold. The correspondence is good although the exact system 
is somewhat less well damped, thus verifying that the effects of the 
zero order hold for a system with the bandwidth of l/4 were ade- 
quately modeled by the first order model of Appendix G. 

E-4 . Comparison With Aided Tracking 

1. Aided tracker description . The application of aided tracking 
to a system similar to the one discussed in this chapter is described 
by Fitts [FI-l] . In this application, the target motion information is 
processed by a digital Kalman filter in which the target is modeled 
as moving in a straight line in inertial coordinates. The output of 
the Kalman filter is a feedforward signal that is added to the feedback 
controller signal. Several versions of this method are described and 
results are presented for tax-gets having both modeled and nonmodeled 
motion. These results are compared with those of a conventional feed- 
back system and found to be much superior. The type of feedback system 
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that is used for comparison is not specified but presumably it is a 
classical-type design since the acceleration error coefficients given 
represent the practical limits that are achievable with such a classical 
design. 


In this Section the aided tracker is compared to the controller 
described in Section D. The application examples by Fitts refer to a sys- 
tem in which the torquer has a natural frequency that is lower than that of 
the system described in this Chapter. However, from the description 
of the aided tracking method, it seems that the results would be no 
different if this method were applied to the system of this Chapter. 

The tracking of a target having a modeled motion is examined 

by Fitts by considering a straight line flyby with a peak angular accel- 

2 

eration of 0.5 rad/sec . Th% results are presented in the form of a rms 
tracking error defined for this particular profile as 


a 


az 



(5.19) 


where is the tracking error. For this target motion, the aided 

tracker has no constant tracking error. The error is determined by 

the detector noise only. This noise is assumed to be o = 5 jlrad. 

d 

For unmodeled target motion, the results are given in the form 
of frequency response curves. Zero detector noise is assumed. 

2. Comparison of feedback controller and aided tracker . Since 
the controller designed above is noise limited, a meaningful comparison 
of this controller with the aided tracker can only be made if the same 
detector noises are assumed for both designs. The comparison of the 
modeled target tracking will therefore be made for two noise levels: 

( a > °d = M-rad, as assumed in Section V-D; (b) g, - 5 [Xrad, as 
assumed in the paper. 

3. — 100 |-lrad . The influence of different detector noise 
levels is examined by Fitts but only noise levels up to 50 |~lrad are 
considered. Extrapolating these results to 100 |lrad, a rms tracking 
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noise of about 70 (trad is obtained. 

For the SEF controller, the rras tracking error is caused by two 
effects: (a) output noise, (b) acceleration error. The numerical 

values were given in Table V-8. Using these values in Eq. (5.19), the 
rms tracking error is obtained as 

0 =75 (JLrad 

az 


4 - QA = 5 [irad . The rms error for the aided tracker is 


°AT = 3,8 tXrad * 


With this detector noise level, the noise limit on our controller is 
removed. The system bandwidth is now limited by the sampling rate. 

The effect of removing the noise limit is considered below for systems 
using DR 2 and D type estimators (see Table V-l). For the DR 2 
system, a high gain controller is used (Design 2) and higher gains 
are selected for estimator No. 1 according to the acceleration error 
equation of App. K. The resulting eigenvalues are shown in Fig. 

V_23 ' 1 Jw (sec -1 ) 


• controller eigenvalues 
O estimator eigenvalues 




h 600 


500 


h 400 


L 300 


200 


b 100 


• a (sec ■*■) 


-600 -500 -400 -300 -200 -100 

FIG. V-23 EIGENVALUES OF HIGH BANDWIDTH SYSTEM 
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The correspondence of the time response of this system with 
that of the exact simulation was found to be satisfactory. Its acceler- 
ation error is 70 |lrad/(rad/sec\ The rms tracking error as defined in 
Eq. (5.19) is 


0 p C = 25 Ur ad. 

The system using the D-type estimator has zero acceleration error (V-D-2), 
It was eliminated for the high noise case because it had unsatisfactory 
disturbance rejection. With the low detector noise, a higher estimator 
bandwidth is obtained and the disturbance rejection is therefore im- 
proved. The maximum deflection for a step disturbance is 85 (irad/lr ad/sec 2 ). 
This value is still higher than that obtainable when the gyro is used 
as an additional measurement (see Table V-8). 

The frequency response curves for unmoceled target motion as 
given by Fitts are reproduced in Fig. V-24. The two figures represent 
two different algorithms. The response of the sampling rate limited 
DR 2 system described above is overlaid on those curves. 

Under the assumption of zero detector noise, the DR 2 system 
is seen to have lower error for frequencies higher than 1 cycle in all 
the cases represented in the paper. When nonzero detector noise is 
considered, the aided tracker has no advantage for even lower frequencies. 

5. Evaluation . For nonmodeled targets, the aided tracker is seen 
to have an advantage over the high bandwidth feedback controller only 
for low frequencies and if a low detector noise is assumed. 

For modeled target motion, comparable errors are obtained in 
both systems when a high detector noise is assumed. For a low detector 
noise, where the system is sampling-rate limited, the aided tracker has 
a lower error unless some disturbance rejection stiffness is sacrificed. 

For this case the aided tracker therefore has an advantage. 

In summary, the high bandwidth controller that can be obtained by using 
full state feedback can match the performance of the aided tracker 
in most cases. The aided tracker has an advantage only in the case of 
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FIG. V-24 FREQUENCY RESPONSE OF AIDED TRACKER [FI-l] 


a target that can be modelled and low detector noise. This advantage 
seems to be outweighed by the added cost of the computing capability and 
by the complexities of the aided tracker. If it is also taken into con- 
sideration that most tracked targets will not follow the model, the 
state feedback approach to the tracking problem seems to give a better 
solution , 
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F. SUMMARY 


The design of a controller for the inner azimuth gimbal of a 
tracking telescope was presented. The controller was designed using 
state estimate feedback. 

To obtain a controller for the entire system, additional controllers 
for the elevation and outer azimuth gimbals are required. The eleva- 
tion controller may be identical to the inner azimuth controller except 
for some minor changes in the numerical values of tho parameters due 
to the different moments of inertia. The outer azimuth controller is 
straightforward and may be designed by classical techniques. 

The special feature of this problem is that the control is applied 
through an elastic element (the hydraulic torquer), which makes it 
difficult to obtain adequate damping at high bandwidths. Still, high 
bandwidth is required in order to obtain a low acceleration error coef- 
ficient. The ability of full state feedback to place the closed loop 
eigenvalues arbitrarily enables the system to obtain the required band- 
width with sufficient damping. The bandwidth and the acceleration error 
coefficient are limited either by the allowable controller noises or by 
the measurement sampling rate, depending on the measurement noise 
levels. These are common limits for controllers designed by state 
feedback. Stability considerations do not limit the bandwidth of such 
controllers as they do for many classical designs. 

The principal design parameters for state feedback controllers 
are the state and control weights. However, important performance cri- 
teria, such as the acceleration error coefficient, cannot be conveniently 
expressed in terms of these weights and only general trends can be es- 
tablished. Several designs are therefore made and their performances 
compared , 

Two measurements are available and it is possible to combine them 
in various ways in the estimator. The comparison between these designs 
is made according to the functional design Criteria. Two measurement 
noise levels were considered. At the higher noise level, the system is 
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noise limited. A detailed comparison between several designs was made 
at this noise level. At the lower noise level the system performance 
is limited by the detector sampling rate. The influence of removing 
the noise limit was considered for two of the above designs. 

The principal conclusions that were reached from the evaluation 
of the various designs are as follows. 

1. The use of the gyroscope as a measurement, in addition to the 
target detector, improved the disturbance rejection properties of the 
system. However, it causes the system to be changed from Type 2 to 

Type 1. 

2. At the high noise level the gyro rate mode is preferable xo 
the rate integrating mode. 

3. The ultimately selected system for this noise level is a 
system using a reduced order estimator. Its performance was made accept- 
able by the use of the sensitivity reduction method as a pole placement 

device. 

4. Removing the noise limit on the bandwidth results in a decrease 
by a factor of 3 in the acceleration error of a system using both the 
gyro and the detector measurement. In a system using the detector measure- 
ment only (which has zero acceleration error), the step disturbance re- 
jection is improved by a factor of 3. 

5. The representation of the sampler and hold as a linear net- 
work proves to be a good approximation for the bandwidths considered, 

as shown by the closeness of the time responses of the real and the linear- 
ized systems. 

6. The use of an aided tracker instead of a high bandwidth feed- 
back controller has an advantage for a modeled target and low detector 
noise only. 

/tty In conclusion, the use of state feedback for this system results in 
a controller that has good transient response and the tracking capabili- 
ties which are comparable to those of a much more complicated aided tracker 
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VI. CONCLUSIONS 


In this thesis, the sensitivity of state feedback controllers (SEF) 
to parameter variations was investigated; a method was developed to 
reduce sensitivity, and it was applied to the Stanford Relativity Satel- 
lite. In addition, the design of a tracking telescope, partially 
aided by the sensitivity reduction method, was studied. 

The principal conclusions in the area of parameter sensitivty 
are as follows, 

1, The use of state estimate feedback, instead of state feedback, 
in systems in which there is a lightly damped elastic element between 
the control and the measurement, results in a consxderable increase in 
the sensitivity to variations of the spring stiffness. Systems of this 
type are common and include structural resonance, . systems with hydraulic 
actuators and others. Their sensitivity is caused by the multiple points 
of unity gain of their open loop transfer function and can be characterized 
by means of the frequency margin of stability. If the elastic element 
is sufficiently damped so that multiple points of unity gain are avoided, 
the sensitivity is considerably reduced. 

2. The sensitivity reduction method that was developed in this 
thesis is an effective design tool, as demonstrated by the results of 
its application to the low order approximation and the full model of 

the Stanford Relativity Satellite. The range of stability was typically in- 
creased by a factor of 2 to 3 . For originally optimal systems, the out- 
put error and the control effort were thereby increased by a factor of *=» 2. 

3. An efficient computational algorithm of the sensitivity reduc- 
tion method is essential for this method to be applicable to fairly large 
systems. Considerable effort was therefore expended in the development 
of such an algorithm. The resulting computer program has reasonably 
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low computation times. The sensitivity reduction of a 12th order system 
requires four minutes approximately on an IBM 360/67 computer. 

4. Sensitivity criteria may be a useful tool for the design of 
estimators. This is demonstrated by the use of the sensitivity reduc- 
tion method for the design of a reduced order estimator for the track- 
ing telescope. 

The following conclusions were derived from the design of the con- 
troller for a large tracking telescope. 

1. The use of a state estimate feedback controller for this sys- 
tem gives it a tracking capability that is comparable to that obtained 
by using a much more complicated aided tracker. 

2. Since arbitrary pole placement is made pooaible by the use of 
full state feedback, there is no stability limit to the system bandwidth. 
The limit is instead determined either by the allowable control level 

or by the target detector sampling rate. Of these two limits, the one 
that is effective depends on the level of the noise inputs. 

3. Two measurements are available and different estimator designs 
are possible using various combinations of these measurements. The 
relative merits of these estimators cannot be determined a priori and 
they are therefore compared according to the performance specifications. 
For noise limited systems, the optimum balance is sought between low 
controller noise on one hand, and low tracking error and good disturbance 
rejection on the other hand, with acceptable parameter sensitivity as 

an additional criterion. 

A reduced order estimator using as measurements, thetarget de- 
tector and the gyro in its caged mode was finally selected for this 
case. 

For the sample rate limited system, the estimator selection is 
more arbitrary. 

4. The use of a saturating controller, together with integral 
control in an inertia-type plant (l/s ), causes an unavoidable large 
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signal instability. This instability occurs at larger signals if the 
input to the integral control integrator is limited. 

5. The modeling of the sampler and zero order hold by a linear 
network is justified in the considered range of bandwidths (up to 
about 40% of the sampling rate). 

In summary, state estimate feedback is found to be an attractive 
design technique for high performance controllers. Its practical 
applicability is enhanced by the use of the sensitivity reduction method. 
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APPENDIX A 


GAIN MARGIN OF STATE FEEDBACK CONTROLLERS 


Consider a single input system 

x - Fx + Cg Gu 
u = -Cx 

where C is a full state feedback gain matrix, and is a variable 

scalar parameter, the nominal value of which is unity. The system will 

have infinite gain margin to changes in c if the roots of the numerator 

-1 

of C(sl - F) G are in the left half plane. This is also true if 
state estimate feedback instead of state feedback is used. 

Proof 

1. State Feedback 

The closed loop characteristic equation of the system is 

det( si - F + C q GC) = det( si - F) det[ I + c Q (sI - F)‘ 1 GC] 

. . = det( si - F) [l + c Q C(sI - F)'^] : 

= det( si - F) + c Q C adj( si - f)g . 

As c increases , m roots of the system tend towards the m roots 

of C Ad j (si - F)G and the others depart in n - m asymptotic direc- 

tions. It will now be shown that 

m — ; n- - 'T::Y 

and that therefore only one root departs on an asymptote, along the 
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The diagonal elements of adj(sl - F) are monic polynomials 
of degree n - 1 in s. This is so because the minors of these element 
are the determinants of (n-1) x (n-1) matrices with s in all their 
diagonal elements. Therefore f even if G has only one nonzero element 
the vector Adj(sl - F)G will have at least one element of order 
n-1 and since C has no zero elements, C adj(sl - F)G is a 
polynomial of degree n-1. If the roots of C adj(sl - F)G are in 
the left half plane, all the roots will therefore remain stable as c^ 
increases. QED. 

2, State Estimate Feedback 

In this case the closed loop characteristic equation is 

det(sl - F + c qGC) det(sl - F + KH) 

= [det(sl - F) + c Q C ad j (si - F)G] det(sl - F + KH). 

The estimator roots are therefore unaffected by changes in c^ and the 
proof remains valid. 



APPENDIX B 


PROGRAM PAROPT 


1 . General 

This Appendix contains a partial listing of the program PAROPT and its 
operating instructions. The listing includes all the subroutines that were 
developed for the program. For the subroutines that were obtained from 
the Program Library of the Computer Science Department at Stanford University, 
only the description is given in the listing. Subroutines from the IBM 
Library (GMPROD and GMTRAN) are not listed. Sufficient comment lines have 
been introduced to make the program self explanatory. 

2. Operating Instructions 

a. Range of the program 

Order of system: 12 

Number of variable parameters: 5 

Number of free parameters: 30 

Number of controls: 3 

bo Input. The program is stored in compiled form. The input to the 
program consists of two parts: (l) subroutine SETUP, (2) Data 

and options . 

(l) subroutine SETUP . This subroutine defines the relationship 
between the parameters (fixed, free, and variable), the 
coefficients of the system matrices, and the control gains. 

It has to be written for each system that is being optimized. 

The input to this subroutine consists of the current values 
of the free parameter vector. The output consists of the 
F, r, and C matrices (see Eqs. 4.3 and 4.9), The subroutine 
is compiled and adjoined to the compiled program. The 
program can then be run with different data sets. A sample 
subroutine SETUP, with the job control language that is 
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required in order to load it on the IBM 360/67 computer of 
Stanford University is shown in Fig. B-l. M405. PAROPT on 
line 39 is the name of the load module that is created. 

This can be changed to any name selected by the user. 

It is recommended that this subroutine be made as general 
as possible so that different sets of parameters can be 
defined as free or variable without changing the subroutine. 

(2) Data and options . The data cards are: 

Card 1 : System description: N, NS, NZ, NP, NC, NW 

N Order of system 
NS Number of variable parameters 

NZ Number of free parameters 

NP Number of fixed parameters 

NC Number of controls 

NW Number of disturbances 

Format : 5 13 

Card 2 : Options: NFLAG, NPRINT, NCO. All these para- 

meters are either 0 or 1 . 

NFLAG = 0. Only initial conditions are required. 

NPRINT = 0. Only condensed information is printed about 
the PI, gradient, and values of the free parameters at 
each step of the search. 

NCO =0: C is not a function of the free parameters. 

Format: 313 

Card 3 : EPS, TOL, . 

EPS The relative weight of the additional cost. If 
EPS < 10" 6 , 6X is not evaluated and only the 
nominal PI minimization problem is solved. 

TOL A measure of the relative precision of the final 

values of the free parameters. Recommended initial 
value: 10" 4 . 

' DIS Use 10" 4 
Format: . 3E12. 5 
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1. // PARO^Tl .inn 'M405, 701,1. , 3', ’MAOASS' 

2. //STF.P1 EXEC FORTHCl, PARM. FORT= , OPT=2' 

3. //FORT . GYT> I M DP * 

I,. 5UP-ROIJT I NE SETI.^C 7, A, R, C) 

5. IMPLICIT RFAl*8 (A-H,0-7) 

n’ DIMENSION 7.(30), A(12,12),B(12,5),C(3,12) 

l\ ■ COMMON/ PARM/ Q( 12, 12 ) , R(5, 5 ), V( 5, 5 )-, QR( 3 ), P( 30 ), S ( 10 )., EPS 


8. 


no 10 1 = 1,6 

9. 


no 10 . 1 = 1,6 

10. 

10 

A( 1, <J)=0. 

11. 


no 20 1=1,5 

12. 

20 

A( 1, l+l)=l. 

13. 


si*sm 

14. 


A( 3, 1)=-P(1)*S1 

15. 


A ( 3, 2 ) =- P( 2 ) *S1-S1 

16. 


A( 3, 3)=-P(3)*Sl-S(2) 

17. 


HO 30 l =1, 3 

18. 

30 

A( 3, 1 +3 ) = D ( 1 )*S1 

19. 


A( 4, 4 )=-Z( 1) 

20, 


A( 5, 4 ) =-Z ( 2 ) 

21. 


A( 6, 2 ) =7 ( 4 ) -S ( 1 ) 

22. 


A{ 6, 3 )=7( 5)-$ ( 2 ) 

23. 


A(6, 4 )=-Z (3) 

24. 


A(6,5)=-Z(4) 

2 5. 


A(6,6)=-Z(5) 

26. 


DO 40 1=1,6 

27. 


DO 40 J=l, 2 

28. 

40 

B( I , «J ) =0 . 

29. 


B( 3, 1)=S1 

30. 


n(6,l)=Sl 

31. 


00 50 1=1,3 

32. 

50 

B( 1 + 3, 2) = -Z(l ) 

33. 


DO CO 1=1,3 

34. 


C(l, 1 ) = P(I ) 

35. 

60 

C(l, » +3)=-P( 1 ) 

36. 


RETURN 

37. 


EKO 


38. /* 

30. //LICFP. SYSUinn np PSNAME=M405. PAPr>' , T,VOLNME = SER=SYSl6,UTnT=2314, 

40. // 0 1 S P= ( MEl’l, KEE R V, SPACE® (TPK, (2,3,1)) 

41. //IKED. a nn DSNAME=M405. PAROPT, Vni.l.!ME=SF.P=SYS16,U?;iT=2314, 

42. // t)tSP»( OLD, KEEP.) 

43. //I.KE0 . SYS I N DP * 

44. INCLUDE 0(?AR0P). 

45. ENTRY MAIN 

46. NAME PAROP 

47. 


FIG. B-l SAMPLE OF SUBROUTINE SETUP FOR SIXTH ORDER SYSTEM 


QV BOOB aBAUW 



Cards 4 and on. 


a. Fixed, variable, and free parameters in this order. Six 
parameters to a card. Each set starts on a new card. 

b. The state weight matrix. Read by rows, six numbers to 
a row. 

c. Control weight vector. For this program the control 
weight matrix is constrained to diagonal form and there- 
fore only the NC diagonal elements are read. 

d. The covariance matrix of the variable parameters. Read 
by rows. Six numbers to a row. 

e. The covariance matrix of the disturbances. Read by 
rows. Six numbers to a row. 

Format : 6E12.5 

A sample data set with the JCL required to run it on the IBM 
360/67 computer is shown in Fig. B-2, M405. PAROPT on line 

2 is to be changed to the name selected when loading the sub" 
routine SETUP. 

(c) Output. For NFLAG = 0, the printout consists of the initial 
dynamics , weight, and covariance matrices, initial eigen- 
values, and initial nominal and additional PI. For NFLAG = 1, 
the printout contains, in addition, search information, as 
well as final values for the free parameters. 

A typical printout for NPRINT = 1 is shown in Fig. B-3. For 
each cost evaluation, the nominal (J 0 ), additional (J A >, and 
total PI (J Q + e.J A ) are printed out as well as the values of 
the free parameters. The total PI is given as a fraction of 
its initial value. At the end of each iteration, more de- 
tailed information is printed out. 

(3) S uggestions for the use of the program . 

(a) Execute an initial run with NFLAG = 0 in order to find the 
initial relative values of the nominal and additional PI. 

(b) Initially select the weighting coefficient, e, and the para- 
meter covariance matrix so that eJ A > 10 J Q . A design that 
is mainly determined by sensitivity consideration results. If 
its nominal properties are unsatisfactory, execute additional 
runs with lower e. 

(c) Limit the time of the run according to the general guidelines 
given in Section IV— C -4 . Even if this causes the run to stop 
before it is finished, the resulting design may be satisfac- 
tory since toward the end the program may become less efficient 
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1 


1. 

//PAr?0 p T2 -ion 

'111*05,301, 1. ,5', 'HAPASS' 

2. 

// J0BL1 B on 

05KAME=Hft05. PAT50PT, VnLL»1E-SER"'>VS16„l)NlT< 

3*. 

// 

ni99*(PLn,PASS) 

ft. 

//STEP! EXEC 

PGH-PAHOP 

5. 

//CT06F001 on 

5Y50HT-A 

c, ■ 

//FTO5F001 on 

* 

7 * 

C 2 5 3 

1 2 

8. 

0 10 0 


9. 

1. 

l.OOE-ft 1.00E-6 

10. 

. ftft 7 

.601 1.216 

11. 

.ft 

.3 

12. 

.88 

.39 -.068 .6 

13. 

.2 



lft, 

15. 

16. 
17. 


18. 

19. 1. 

20. .16 

21 . 1 , 

22 . /* 


.09 

1 . 


FIG. B-2 SAMPLE OF DATA FOR SIXTH ORDER SYSTEM 


Results of 1 f 
PI evaluation^ 


0.10 25 ft 537 69247660 05 
0. in c ftoi0 t )9Al<>750 02 
-0.3075015172653910 02 


9.1025058769252710 05 
free i 0. 10.556910996 19750 02 
parameters * -0.3970015172653930 cz 
0.102585 B76924506D 05 
0. 1355491099819750 02 
-0.3079015171163810 02 
0.l025e5tl7fc9337830 05 
J 0.1855451 09 9619750 02 

O -P.3079O16172653930 02 


STAI’JS AT 1 TER AT'ItIN H 1 


0 •5474052 762 24 6790 03 

0.1752326903699040 03 
0.2622043650432460 02 
0.547 605? 755 60 50 3D 03 

0.1752326903553430 03 
0.2622043698432460 02 

0.5476052762555480 03 

0.1742326903550930 03 
0. 2622043658432460 02 
0.5476052767849180 03 

0.1752326903550830 03 
0.2622043659922480 02 



0.9495198503304340 00 

0.1546&69666465040 02 

0.9495198502 3R0430 00 

0.1546659667955160 02 

t<. 949 5 19 8 50 85 72 6 70 00 

0. 1546669666465040. 02 


0.9495198513723950 00 

,1546669666465040 02 


\ 0.154 

\ J + 


eJ. (relative) 
A 


CURRENT solution 
0. 1856491 9 QQ 6l975f> 02 
0. 1 75232690355083D 03 
0.1546665666465040 02 
-0.3079015172653930 02 
0.2622043698432460 02 


GRA0IENT 

0. 1C36 202 535033230-02 
0.4331729374325950-02 
-0.3542291931803000-01 
0.6132469581199650-02 
0.4070217069238420-01 


DIRECTION OF SEARCH 
-9.17319531 01371570-3 1 
-0.2336833264425980-0 l 
0.2132676036914590 00 
-0.3903200790824450-01 
0.9217989013351290-01 


APPROXIMATE MINIMUM VALUE OF ft 0.9495198507658860 CO AtFKA: 0. 76743593100 02 

NUMBER OF FUNCTION EVALUATIONS: 16 


FIG. B-3 TYPICAL PRINTOUT OF PROGRAM PAROPT 


LISTING OF PROGRAM PAROPT 


C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


PROGRAM PAROPT 
IMPLICIT P F Al*A I A -H, 0-Z ) 

LOG I C A I P ° I NT, PRTFC 

DIMENSION II 301 

C CM MON /C R I TFR/ ETA, TOL, STfcPMX, OpPS 
COM mom /LOG! CL/ BRIEF, PRINT, UNITL 
EXTERNAL COST ,STAPTR 

510 PQP VA T ( 613 ) 

511 FORMAT ( fiFt 2 .51 
513 FORMAT (AIT) 

M2 FORMAT ( 3F12.5) . , _ n lurTrr , c , , 

3 55 p qp M A T ( • 0 • • 'ORDER CF SYST p M = * , I 2. / VARIABLE PARAMETER- *12,/, 

IFPCf pApAf'gTFRS -* * T? , /, ' FP5=»,F7.3) 

756 F- 0° MA T < • (' • , * I NT T I A L COST ONLY R ECU TPFD* I 
3 C 7 pppyATI ' * , 6 X, • NCC* = ' , 12 ) 

M IS THE ORDER OF THE SYSTEM 
NS IS THE MUMREP OF THE V AR I A 8L F PARAMETERS 
Nl IS THE NUM8 p R OF FREE PARAMETERS 
NR IS THF NUMBER OF FIYED PARAMETERS 
NC IS THE NUMBEP OF CCNTRCLS 
MW IS THE NUMBEP OF 0 IS T URRA NCF S 


READ! 5,510 N, NS, NZ , NP , NC , NW 

IF N FI, AG= 0 ONLY INITIAL CONDITIONS ARE REQUIRED 

jF vpRlNT=0 SEARCH PP INTOUT IS CONDENSFD 

IF NCCI=0 C IS NOT A FUNCTION OF THE FREE PARAMETERS 

RFAD(5,513J NFLAG, NPRINT, NCO ,NXOF 

EPS IS THE PFLATIVE WEIGHTING ON THE ADDITIONAL COST 
TOL IS TH r - RELATIVE PRECISION fl F THE FINAL VALUES OF 

THE EPFF PAPAME T ERS „ 

DIS IS THE RFLATTVE pFRTUFBATION OF THF VARIABLE PARAMETERS 

READ (5,5121 FPS, TOL, DIS 

PAR- FIXED "PARAMETERS . 

SI - VARIABLE PARAMETERS 
l - FREE PARAMETER S 

READ! 5,511) (PAR ( I ) ,T =1 *NP V 
P C AD ( 5 v 5 1 ! 1 ( S I ( 1 1 , I - 1 , N $ ) 

READ! 5 , f l IMZII 1 ,1 -1 , NZ » ' 

0 - THE STATE WEIGHT MATRIX 
QR- THE CONTROL WEIGHT VECTOR 
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P - THF COVARTANCF MATRIX OF THF VAPIAPLF PARAMETERS 
V - THE C'OVAMAMCE v AT0IX OF Tup OISTUPflANCFS 

READ! *,511) 1(01 t, J>,J=1 ,N> ,| = l ,N> 

9 FA D! 5 * 51 \ ) (CRf I ) , 1 = 1 ,NC) 

READ! 5 » 51 1 1 l 1 0 ( I * J ) , J =1 ,NS ) » I =1 ,N S I 
READ! 5. ^111 ( (V! T » J 1 , J = l ,*fW ) , ! = l ,N W I 
WRITE If ,355 )N,N$, NZ.EPS 
WR I TE ( 6 * 3 C 7 ) NCO 

IFINFI AG.FO. 1130 TO 323 
W3ITE ( * ,3561 

323 I F I N?P I NT , EO, 01 GO TO 1 Of; 

PR ! NT = • TRUE » 

BFI FF =. FALSE, 

GO TO J 10 
100 PR INT= , FALSE* 

BRIEF-. TRUE. 

110 CALL DIN I T (COST ,NZ ,Z ) 

IF CNFLAG.FQ.OIGO TO 600 
NL D I MsMZ* (NZ—1 > /2 

CALL MINMZFtNZ, Z, NLDIM, COST, STARTR) 

600 STOP 
FNO 


OMGlNAIi PAGEIS 
OF POOR QUALITY 
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I 


SU'VRniJT IAJ E DINIT J 

THIS SUnoPUTIMp WRITES TH«= INITIAL DATA AMD C ALCUt A TF.S THE | 

INITIAL PI AMD " I r.^MVAL U r - S USING SUBRHUT! N K COST ( 

INPUT TO SUBROUTINE DINIT . i 

NZ : THP NUMBER OF FREE PARAMETERS | 

l : THE FR?F PARAMETER V C CTDP i 

FUN : A SUBVnUTTNF- RETURNING THE VALUE OF VHP PI AT THF | 

POINT DEFINED PY THE VECTOR Z. i 


SUBROUTINE OINI T( FUN, N 7 , Z! 

IMPLICIT REALMS (A-H,0-Z» 

DIMENSION A< 1 2 , 12 ) ,B( 12 ,E ) , 7 ( 3 D ) ,C < 3,1 2 ! 

C n.Y MO N N , NS , N'P , NW , MC , NC 0 
C ON Mr) N / po q t s / C R (l?) f CK!2) 

00M V0N/ PARM / Q( I 2 , I ? ) ,R ( f , «? ) ,V( S , 5 ) ,QP( 3) , PAR (30 ! ,S I ( ID I , EPS , 0 IS 
COMV.rjN/TO/'FPl ,& Pi 
C DM NON /MOPES /MD QF 
CPMMQM/ FV 4 L/NEV AL 

2 n l FQRm a t( i.;ii , //,?X, • the INITIAL DYNAMICS MATRIX.,.'! 

20 R FO»MAT( * 0 ', //, 2 X, • THF DISTURBANCE DISTR I 9 UT ION MATRIX...'! 

21 A F C'R.MAT ( 'O' , / / , 3 0 X , 'INITIAL COST',//, IPX,' BASIC COST'.llX, 'ADDITION 
1 AL COST 12 X, TOTAL COST'! 

202 FORM AT ( 'a*, 4 X, 3 < 1PLP5.15) 1 

20 4 FORMAT ( • 0 • , / / , 3 OX , • IN! T ! A (. PHOTS' , //, 20 X, 'RFAL ' , 24 X , • I MAG .* ! 

216 FOR MAT ( • O' , //, 2 X, ' THF DISTUPBANCF COVARIANCE MATRIX...*! 

2 C 6 FOP MAT ( 2 ( 1 3 X * FI A ,8 I ) 

4 HC FOP V.AT( •:;»,//, ex,' THE PAPAMF*- = c covariance MATRIX... •» 

ICO FORMAT ( * 0 ',//, 2 X, ' THE STATE WEIGHT...'! 

IG 1 FORMAT ( 'O' , / /,2 X, 'THE CONTROL W c ICHT...'l 
CALL SFtup( Z, A ,B, C ) 

WRITE ( 6 , 201 ) 

CALL v A TO U T ( A , 1 2 , 1 2 , N , N I 
WR I TE ( 6 , 100 ) 

CALL MATOUTIQ, 12 , 12 , N, N I 
W C ITE( 6 , 101 I 

W P I TE ( 6 , 2 (' 2 ! ( 03 (I ! , I = 1 , N'C ! 

WRITFf a, 20 B> 

CAL L M A TOUT ( 3 , 12 , 5 , N , NW ! 

WRIT" ( 6 , 216 ! 

CALL MATOUTI V, 5 , 5 ,NW,NW! 

WPT TE ( 6 , AQO I 

CALL M A TOUT (R , 5 , 5 , NS, N$ 1 

neval=o 
w 0 PF =0 

CALL F UN (HZ, Z , F I 
W R I T F ( 6 , 2 1 4 ) \ 

WP 1 TE (6 ,e r ?! FPI ,API, F 
WRITE( 6 ,?C 4 > 

WFITE ( 6 , 206 ! (CR II! ,CI ( I 1 , 1 = 1 , N) 

RETURN 

END 
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SUBROUTINE COST 


THIS SUBROUTINE CALCULATES THE COST USING SUBROTINE OERIV 
FOR THE CALCULATION OF X ANT DX. 

INPUT TO SUBROUTINE COST: 

N Z : THF NUMBER OF F REE PARAMETERS 

Z : THE FR C E PAR AN'FTfR VECTOR 

OUTPUT FROM SUBROUTINE COST 

PI : THE VALUE OF THE PI AT THE POINT DEFINED BY THE 

VECTOR Z. 


SUBROUTINE COSTINZ, 7 , PI) 

IMPLICIT PEAL*8< A-H.O-Z ) 

DIMENSION 7(30) ,XX ( 12 , J 7 ) , 0X( 12 , 1 2 ) , B M 2 , 5 ) ,01 ( 12 , 1 2 1 , 

1 C ( 3 * 1 2 > 

COMMON N f NS » NP » MW * NC » NC 0 

COM MON /PAR M/ Q( 1 ? , 12 ) ,P I 5 , 5 ) , V ( 5 , 5 > , 0R( 3 ) ,P AR ( 30) ,S K 10 ) , EPS, DIS 
COMMON /TQ/FP I ,A Pf 
CO^MON/EVAL /NEVAL 

COMMON/MODES /mhoe 

C0MM0N/ST0PE/GXU2 ,12 ) , GCX ( 1 2 , 1 2 ) , GY ( 12 , 12 , 5 ) , GYY (1 2, 1 2 , 5 ) 

220 FORMAT ( »0 ’,//, 2 X, > THE STATE COVARIANCE MATRIX*) 

CALL DERI V(XX,DX,Z ,C, CfcSC ) 

F P I =0 • 

API =0. 

I F( NEVAL# GT« Oo AND* NCO» E0» 0 ) GO TO 41 
DO 54 1=1, N 
DO 54 J=1 , N 
5A Oil I, J )=0. 

DO 55 1=1, N 
DO 55 J = 1 , N 
DO 55 K=l ,NC 

55 Ql( I, JI=C(K,T)*C( K,J)*OP(K)*Q(I ,J) 

41 IF(NEVAL,GT.O)GO TO 150 
F 1 = 1 . 1 

WR IT E ( 6, 220 ) 

CALL MATOJT( XX, 12,12, N,N) 

150 DO 50 J =1 , N . 

DO 50 K=1 , N 

FPI=FP?>Qi(K, J ) *XX< J,K> 

IF IB PS.l E • 1 • 0- 6 ) GO TO 50 
AP I =AP I *Q 1 ( K ,U) $ DX ( J , K ) 

50 CONTINUE 
GO TO 60 
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650 FPI=l,D6*rPI 
A P I =0 . 

60 Pr=FPI+EPS*API 

IF(NPVAL.GT.O >G0 TO 68 
FO = P I 
GO TO 65 
68 PI=PT/FO 

I F ( MO OF • F Q« 1 . OR • PI • GT , F l I GO TO 65 
FI = p I 

DO 66 !=l t N 
DO 66 J=I ,N 
GX( I,J > = XX ( I ,J) 

GDX{ I, J)=QX< !, J) 

GX( J , I ) =GX( I , J» 

66 GDX ( Jt I ) = GDX ( I , J ) 

DO 67 I = 1 » N 

DO 67 J = 1 , N 
•DO 67 K= 1 ,N$ | 

67 GYV( T » J , K > = G Y{ I,J,K> 

65 WRI TE (6,208) FPT ,API,P 1 
208 FORMAT ( 1( 6-X , 025 » 1 5 ) ) 

WRITE(6,?12) (ZC I ) , 1=1, NZ) 

212 FORMAT ( "»( 6X ,D?5. 15 ) 1 
NEV AL=NEVAL+1 
RETURN 
END 
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SU9 ROUT I N F PER I V 

THIS StJ'.'SfTUT INE DETERMINES THE VALUES D F X AMP DX FOR 4 
GIVEN 7 VECTOR. 

INPUT TO SUBPOJTINE OEPIV: 

Z : thf free parameter VECTOR, 

OUTPUT FROM SUBROUTINE DERIVt 


XX : THE NOMINAL COVARIANCE MATRIX. 

DX : THE ADDITIONAL COVARIANCE MATRIX. 
C : T HE CONTROL GAIN MATRIX. 


SUBROUTINE DE° I VI XX, D X, Z ,C , *) 

IMPLICIT REALMS (A-H.n-Z ) 

DIMENSION XXI 12 ,12 ),DXI12,12> ,RU 2,5>,BV( IZ,5> ,BT{*,12> ♦ A(12, 12), 

1 AH12, 12) ,B1I12 ,5 ) ,C( 3, 12 ) ,C1 (3,1?>,0A( 12,12) ,DR(1 2,5) ,0AT(l2 ,12) , 
MJL< 73 ,70) ,1 PSI7 8I , C3M2 , 12 ) ,DAX 1 12,1? ) , ULU7fl ,78 ) , I PS I I 78 ), 

3QE17P) ,COV( 75), Y( 12, 12) ,71 301 ,D< 12), INT( 1 2) ,GA< 12,12) r 
AG AX ( 12,12) 

COMMON N, NS , NP , NO, NC , NC 0 

COMMON/ parm/ 0( 12, 12 ) ,R ( 5 , 5 ) , V ( 6 , 5 ) , QR< 3) , P AR ( 30 ) , S 1 1 1C) , EPS, 01$ 
COMMON/ROOTS/CWRI 1 2) . CW I ( 12 ) 

COMMON /MODES /MODE 

COMMON/STORE /C.X (12,1?) ,GOX< 12, 12), GY (12, 12 ,5) ,GYYI 12,12,5) 

CALL Sf TUPI Z, A, B, C ) 

IF ( MODE. 60.1 )G0 TO 201 

MODE =?- THE INITIAL VALUE Op THE PI OP THE PI ALONG A LINEAR 
SEARCH point is EVAIUATED. 

MODE =1-TH C PI IS EVALUATED FOR GRADIENT DETERMINATION. 

PER TUB AT I ON EQUATIONS APE USED. 

DO 110 1=1, N 
DO 110 J=1,N 

G A ( I , J ) = A ( I , J) 

110 Al l I , J 1 = A ( I , J J 

EIGENVALUE DETERMINATION. USFO TO DETERMINE STABILITY OF 
SELECTED POINT. 

CALL BALANC (12, N , A 1, LO, I H l , D ) 

CALL FLMHES 1 12, N, LO, I HI, Al , I NT ) 

CALL HQR (12, N, LO, I H I » Al, CWR, CWI, I ERR ) 

DO ?O0 i=1 t n 

IFICWR (I ) .GE.O. ) GO TO 720 


2CC CONTINUE 
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C EVALUATION OF THF RIGHT HAN 0 SIDE OF THE EQUATION 

C 

C (U A*XX+XX*A 1 =-B*V*B 1 > 

c 

C THE RIGHT HAND SIDE IS TEMPORARILY STORED IN XX. 

C , 

20 1 CALL GMPono ( B,V , 8V,N,NG,NG, 12 ,5 > 

CALL GMTRAM<R t 8T,N,NG,12,5) 

CALL GMPR OD ( BV, BT» XX, N, NG *N ,1 2 * 5) 

I F ( MODE. EQ .0 ) GO TO 202 
C 

C EVALUATION OF THE ADDED TERMS OF THE RIGHT HAND SIDE OF 

C THE PERTURBATION EQUATION 

C 

C ( 2 1 A*XX + XX*A* = -B*V*B , -A*GX-GX*A». 

C 

C WHERE GX IS THE COVAP IANCF MATRIX AT THE NOMINAL POINT AND 

C XX IS THE PERTURBATION OF THF COVARIANCE MATRIX DUE TO THF 

C PERTURBATION OF Z. 

C 

CALL GMPRODI A,GX,GAX,N,N,N, 12,12> 

DO 203 I =1 , N 
DO 203 J = l ,N 

XX ( I , J)=XXt I , J I +GAXII, J)>GAX<J , I) 

203 XXI J. I )=XXI I , J) 

CALL SCOVIGA , XX ,UL,QE,COV, IPS, 2) 

DO 2 DA l-l , N 
DO 20A J = I *N 
XX( I ,J)=XXU * J ) +GX (I » J I 

204 X X ( J * I > = X X 1 1 , J ) 

GO TO 205 

202 CALL SCOVIGA, XX, UL,QE,CCV,!PS,1» 

205 DO 5 1=1, N 
00 5 J=1,N 

5 0XII,J>=0. 

IF! EPS.LT.l. 9-6 ) GO TO 100 
C 

C EVALUATION OF DA/DS 1 1 1 AND DBVDS I1 1. 

C 

DO 90 M=t ,NS 
OS I =D I S*S I ( M ) 

S ID=S I ( Ml ' 
si(M) = snM)«-Dsr 
CALL SETiJPIZ,Ai ,B1 fCU 
■ SUM I -S I D W'-'V 1 
DO 30 1=1, N 
DO ID J =1 , N 

10 0 A ( r , J > - < A 1 < I , J > — A ( T 

or 2 d k =i , no 

20 D B I I , K > = I B 1 1 I , K I - B ( I , K ) » / D S I 
30 CONTINUE 
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EVALUATION PF THE RIGHT HANO SIDE QF THE SYMMETRIC EQUATION 

<31 A ’f'Y + y * A ' = -B # V* ■( DB/DSI l> )- XX* (DA/0 S< l>» 

CALI GMTR AN ( DB » BT ( N, NG, 1 2,5) 

CALL GMPRDCM RV, BT , A1 ,N, NG, N,1 2, 5» 

CALL GMTRAN(PA,DAT t N,N, 12,12) 

CALL GMpO'ODtXX* DAT , DA X, N, N * N, 12*121 
DD 40 I =t t N 
DD 40 J = I » N 

C Q( I . J) =DAX ( I , J ) +OAX< J, I ) + Al( I , J > 4 AT < J, I) 

-40 CQ( J, I ) =CQ< I , J) 

IF < MODE. EQ.O I GO TO 207 

EVALUATION OF THE ADDITIONAL TERMS OF THE SYMMFTR IC 
PERTURBATION F QU AT I ON 

< 4) A*DY4DY*A 1 =-A*GY-GY*A *-B* V* ( OR/DS (I ) >-GX* ( DA/OSI I ) ) 

WHERE GY IS THE VALUE OF Y AT THE NOMINAL POINT AND 
DY IS THE PERTURBATION IN Y DUE TO THF PERTURBATION IN l. 



DO 

203 

1 = 1, N 


DO 

208 

J = 1,N 

20 B 

GA X ( I, 

J>=0. 


DO 

209 

1 = 1, N 


DO 

?C9 

J = 1 ,N 


DO 

209 

K = 1 , V 


200 GA X ( I , J ) = A { I ,K)*GYY( K» J,M)+GYY ( I , K , M ) * A ( J, K > +GAX C l , J ) 

DO 210 1=1, N 
DO 210 J = 1 » N 

210 COU , J > = CQ ( I . J ) +GAX( I, J ) *0 AX (J » I > 

SOLUTION OF THE SYMMETRIC EQUATION (3), 

207 CALL SCOVI A,CO,UL,QE,COV, T PS , 2 I 
DO 50 I -1 »N 
DO 50 0=1 ,N 
50 YU , J ) =C0< I , J) 

EVALUATION OF THE RIGHT SIDE OF THE ANTISYMMETRIC EQUATION* -II. 
Nl=N-1 

DO 60 1=1 ,Ni 
11=141 

DO 60 J = I1 »N 

CO I I , J ) = 0 AX ( I , J 14 All I , J ) -HA Xt J * I ) — A 1 ( J , I ) 

I p ( MODE. F.QoO >00 TO 6 0 

EVA l NATION OF THE ADDITIONAL TERMS OF THE ANTISYMMETRIC 
EQUATION! A). 



ooonoooo onnnoonoon n o o r> o o 


CQ( I ,J)=CQ(I ,J)*GAX( I,J)-GAX( J, I) 

60 CONTINUE 

I F ( NODF • EQ ,0 ) GO TO 211 

SOLUTION OF THE ANT IS YMMFTR I C EQUATION (41. 

CALL SCOVA( A »CQ ,UL1, QE ,CCV, I PS 1,21 
GO TO 212 

SOLUTION OF The ANTISYMMETRIC EQUATION! 3>. 

211 CALL SC OVA ( A ,CQ »UL 1 , QE » COV » IPS l »M I 

212 DO 7C 1=1 ,N 
DO 70 ,N 

70 Y ( I » J ) = ( C. Q( I » JT +Y ( I , J H /2 • 

I F < MDDE.FQ.i ) GO TO 213 
DO 214 1 = 1, N 
DO 214 J=l,N 

214 GY ( I * J *M) =Y( I » Jf 
GO TO 215 

213 00 216 I'=1,N 
00 214 J - 1 * M 

216 Y(I,J>=Y(I,JM-GYY( I,J,M) 

EVALUATION OF THE RIGHT HAND SIDE OF THE FQUATION 
(51 A*OX + OX*A * = -SI GMA t DAT (M d A1 ( M) *»{ «, Ml) , M=1 ,NS 

WHEPE 

DAT (M)=(DA/PS(M) ) *Y 

A1(M) = (D3/DS( M ) ) * V*( DB/D S ( M ) )• 

NS IS THE NUMBER OF VARIABLE PARAMETERS. 

215 CALL G” PR OD ( DA, Y,0 AT, N,N,N, 12 * 1 2 > 

CALL SYMPRD(DB, V, A1 , NG, N, 1 2 ,5 » 

DO 80 1=1, N 
00 80 J = ! , N 

DAT ( I ,J ) = OAT ( I, J ) ♦ DAT (J »I I 
80 DAT ( J , I ) = C'AT ( I, Jl 
DO 85 1=1, N 
DO 85 J=1 ,N 

85 DXI I ,J> =( DAT! I, JJ+Al ( I , J.1 ) *p ( M, M I ♦ DX ( I , J I 

00 CONTINUE ■ - 

IF ( MODE . FQ. r ) GO TD 221 

EVALUATION OF THE ADDITIONAL TERMS OF THE RIGHT HAND SIDE 
OP THE P'ER HP RATION EQUATION 

( 5 ) A *POX + PDX fr'A ' =-S I G MAC 0 AT f M )> A 1 ( M V*.fl ( M , M > J - A*OOX-GOX* A* 

. ' M=1 »NS ' 

WUFRE GpX IS THE VALUE OF OX AT THF NOMINAL POINT, 


n r> o 


call r.MpP> n r.( a , rnx * oa x » m »n»n , 1 ? » 12 » 
no ?i7 1 - 1 » n 
0"! 217 J = J i N 

?i 7 nxn ,j ) = ox n , j) +r,AX( i,j» fGAXu.n 

SOLUTION OF EQUATIONS ( 03 ( f?) * 

221 CALI. SCnVI A,nx.ULfQE .COV.iPSi2 > 

IF ( MODF. EQ,0 I GO TO 10? 

21 B DO 2»0 I— 1 *M 
no 22? J = l tM 

220 DX( I , J ) = OX(l # J ) *GOX( I » J > 

GO TO 100 

720 OF TURN 1 

100 SFTIRN 
END 
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SUBROUTINE SCOV 


THIS SUBROUTINE SOLVES THE EQUATION 


A *X+X*4 • =G 
WITH SYMMETRIC G. 

THE EQUATION IS TRANSFORMED INTO THE FORM 
V *Z =B 
WHERE 

7 IS THE ( N+'i) *N VICTOR OF THE EIF.MF.MTS OF X, 

R IS T HE CN+5)*N VECTOR OF THE ELEMENTS OF G, 

F IS OBTAINED FROM A PY THE TRANSFORMATION OF THIS 
SUB ROUT INE. 

INPUT TO SUBROUTINE SCOV: 

Q : THE MATRIX G 

A VCN : THE MATRIX A 

MODE J MODF = I-v IS COMPUTED! EVAI UATION OF XX) 

MQOFM-THE STORED V IS USED! EVA LUAT ION OF Y AMD OX) 

OUTPUT FROM SUBROUTINE SCOV: 

UL : THE LU TRANSFORMATION OF THE V MATRIX 
QE : THE VEC TOR 0 

COV : THE VECTOR Z 

Q : THE M AT P I X X.NOTE THAT THE MATRIX G IS NOT PRESERVED. 


I 

1 

I 



SUBRQUT INF SCOV! AVCN, 0, UL , OE, COV, IPS, MODE) 

I MPLI C I T R E A L*R ( A - H, 0-7 ) 

01 MENS ION AVCN! 1 2 , 1 2 ) , Q ! 1 2 ♦ 12 ) , CO V! 78 ) , QF ( 78) , 
1UL ( 7B » 7 R ) ,L ( IP, 1 2 ), IPS! 78 ) 

C DM MOM / VM A T / V! T B . 7 B) 

COMMON N.NS.NP, NGf NC* NCO 
M=N*!N«-l)/2 

k = o 

•DO 2.1 I = 1 , M 
■ DP ?Q J = !,N 
: K \ = K+l 

QR<K) =-Q ! I » J) 

L ( I » J ) = K 
20 L! J,I ) = K 

l T ( «0 OE ,= 0.2 )G0 TP 90 
. DO ?-v t = 1 9 m 
00 30 J = 1 . M 
30 V ( I »J ) =0.00 
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on '*0 I = l»N 
n,o 6 0 j = l . 'i 
on 

'.c VILI!,KI,t.(JjK!l = AVCM U , J )♦ V( L( I ,K ) , L < J , K )) 

on 50 i = l * n 

on J = 1,1-* 

50 V<L(I t I f, I» f J> 

qp C At L 1 1 N’SY? ( mot f* M ,V, 7R t CE ,rnv» UL i I PS VO IQ ITS* f. 1M ♦ f.l *2 , £! 43 > 

K = 0 

00 80 I = l . N 
On 80 J = I , N 
K = K«-i 

o 1 1 , j > <nv< k> 

80 3 ( J*» I ) = Q( !» J) 

GO TO 155 

i4i wr net f-, i5i) 

GO TO 1-55 
16? IT= ( c» l 52 ) 

GO TO 155 
IV 3 w P I T £ ( r » 1 *3) 

GG TO 155 

l r l FQ R M A T ( * * ** * M A T 3 IX WITH A POW OF ALL Z C PQ ELEMENTS FOUND* ) 

15 2 format ( PIVOT ELEMENT FOUND*) 

153 FO° mat ( *>l , ***MAX I MUM NUMBER OF T TEF ATI ON $ REACHED IN I M PROVE*> 
155 CONTINUE 
RETURN 
END 
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subroutine sccva 

THIS SUBKCUT I NE SELVES THE EQUATION 
A*X+X*A' =G 

hi TH ANTISYMMETRIC G. 

THE EOLATION IS TRANSFORMED INTC THE FORM 

V* Z=B 

hHERE 

Z IS THE ( N-l ) *N VECTOR OF THE ELEMENTS OF X, 

B IS THE ( N— 1 1 *N VECTOR OF THE ELEMENTS OF G, 

F IS OBTAINED FROM A BY THE TRANSFORMATION OF ThIS 
SUBROUTINE. 

INPUT TO SUBROUTINE SCCVA: 

CQ : THE MATRIX G 

A : THE MATRIX A 

M : M=i-V IS COMPUTED (FIRST VARIABLE PARAMETER) 

M> l-T.HE STORED V IS USED (OTHER VARIABLE PARAMETERS) 

OUTPUT FROM SUBROUTINE SCOVA: 

UL : THE LU TRANSFORMATION CF THE V MATRIX 
QE : THE VECTOR B 

COV : THE VECTOR Z 

CQ : THE MATRIX X.NOTE THAT THE MATRIX G-IS NOT PRESERVED, 


SUBPOUT INE SCPV M A ,CQ , UL , Q c ,C OV , I P S, M ) 

IMPLICIT RFAL*B( A-H, D-Z ) 

DIMENSION A (12# 1 2 1 » C 0 1 12 * t 2 ) , UL ( 7 8 ,78 ) , IPS17RV , 
1 QF( TP) , C P V ( 7 p ) 

COMMON N , NS , MP, MG , NC , MC 0 
COMMON AVM AT/V ( 7 B ,7.6 ) 

MN = N’< t ( N-l ) /?. 

DO 70 ! =1 , MN 
DO 70 J = 1 »MM 
to Vl I *J v-o. 

: K-0 
N-l -N-l 

DO PC 1=1 i N1 
11=1+1 

OP 80 J =! 1 , N 
K=K + 1 

80 QE (K> =-C0 ( I ,.n 
M OD E = 1 

I F ( M • G T • 1 ) GO TO ICO 
NT = N— ]. 

00 1 5 1 =1 , N 1 
00 10 .1 = 1 , N 1 
.1 r v ( I , J ) = A (1 vl ,.i + 1) 

15 V( I i I > "VI I , I 1+1 (1,1V 
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DO 60 K =2 * N! 

J1=(K-1>M 2*N-< >/? 

K1=K-1 

NK=N-K 

L=o c 

DO 30 I =1 » K1 
DO 20 J -I »MK 
VtKl+L * Jl ♦ J 1 =-A( I t J+K ) 

VIKHL* J. JU J)=A( I tK> 

V(J1+J»K1+L)=-A< J+K, ! ) 

VCJl + J»Kl+L*-J)=a(K ,!> 

20 CONTINU p 
L=L+N-!-l 
30 CONTI NUE 

DO 50 I=1,NK 
DO &0 .1=1 1 NK 

40 VIJIM . Jl+J) = M I«-K,J+K) 

50 V<J 1+ l t JIM > =4< K,K> + V(J1M . JIM > 

60 CONTINUE 
GO TO 110 
ICO MODE=2 

110 CALL LI MS V 2 < MOD E , MN, V . 7 8, QE »CCVf UL t IPSt DI GITS ♦ £141 » £1 42 » f, 1,4? ) 
K=0 

Nl=M-l 

DO 90 I =1 « N 1 
11=1+1 

DO 90 J=I1'.N 
K=K+l 

CO( I * J ) =C OVIK) 

CQ ( I, I I =0 

90 CQ< J» I I =-CQ( I ,J ) 

C Q I N » M = 0 • 

GO TO I 55 

141 WPI TF (6 , 1 51 > 

GO TO 155 

142 WRITE <6,1 52 I 
GO TO 155 

143 WR I T E ( 6 » 1 53 J 
GO TO 15 c 

151 F DRMATJ if P. I X WITH A ROW OF ALL ZERO ELEMENTS rOUND* J 

152 FORMAT! «****ZF*n PIVOT FI E MEN T FOUND' ) 

153 FORMAT ( ' * * ♦ * M A X I MUM NUMBER 0 P ITERATIONS REACHED IN IMPROVE* I 
155 CONTI NUF 

RETURN 

end 
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S UR ROUT IMF MATOUT 

THIS SU3RCUT INE WPITFS THE MATRIX A 
INPUT TO SUBROJTINF MATOUT: 

A : THE MAT R IX A I l*J) 

MO! : THE NUMRFp OF ROWS IN T HF STORAGE ARRAY OF A. 

N02 : THE MUM RE P OF COLUMNS IN THE STORAGE ARRAY OF A. 


SUBROUTINE MATOUT! A ,NO! ,N0 2,I ,J> 

0 I MENS ION A! NO I , NO 2 ) 

OnUBL F PRECISION A 
! FORMAT ( r H9 ®OW f 7! PX,AKOL . ,13,1X1) 
2 FORMAT (IA,^X, IP 701 (S.B ) 

DO 10 JJ=l,J,7 
J7=JJ*6 

IF(J7.GT.JI J7= J 

WRI TE I E ,1 ) I JK,JK=JJ,J7) 

00 10 I 1=1, I 

10 WRITE (6, 2 )I I ,(A ( I I ,JK ), JK=JJ»J7) 
RETURN 
END 
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SUB ROUT INF SYMPRD 

THIS SUBROUTINE FINDS THE PRODUCT 
C = A*B *A ' 

INPUT TO SUBROUTINE SYMPRD : 


A 

B 

11 

12 


the matrix A(N*M» 

THE NUMBER O^POWS IN THE STORAGE ARRAY OF A 
;st SJSIE Sf columns IN THE ST ca age of * 


OUTPUT FROM SUBROUTINE SYMPRD: 
C : THE MATRIX C!N*N> 


10 


SUBROUTINE SYMP P D( A»B»C»M»N»I1 * l?) 

I ° L I C I T PFAL^BIA-HrO-Z) 

01 M Ef'S I n N A l 1 1» M ) t B( I 2 » M I »C ( 1 1 * N) »D( \ 2 » 1 2 I 
DO 10 1=1 »N 
DO 10 J-ltN 
D('I »J) =0. 

C( I r J ) = 0. 


on 2 0 1 =1 t N 
DO 20 J = 1 » M 
DO 20 K =1 *M 

20 D( !,J >=M I»K1*3 (K, J>*0( !»J> 

nn ao i =i » n 

DO ao J = I * N 
DO AO K =1 f M 

AO C ( I . J >=D( 1 ,k)*VIJ,K)+ri !,J> 
DO 30 1=2 tN 
T N = I *• 1 

DO 30 J=1 » I N 
3 1 . C (I t J ) ~ 0 ( J ♦ 1 1 
.RETURN s 
END 
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c ------ 1 

c 

SURROUT INE L INS Y2 < MODE, N, A, IOIM, B, X, IU, IPS, DIGIT'S. »*.,*,*) 

C PROGRAMMER: S. R. OUTROW, JR. 

C COMPUTER SCtFMCF OFP A» TMENT , STANFORD UNIVERSITY 

C CONVERTED FOR USE WITH THE WATFOR COMPILER FEDRUAPY 1970 
C MICHAEL MALCCLM, COMPUTER SCIENCE DEPARTMENT, STANFORO 

C 


BASED ON THE SINGLE PRECISION REAL LINEAR SYSTEM 
PACKAGE WRITTEN BY JOHN W. WELCH OF SIAC. WHICH IS 
IN TURN, AN ADAPTATION OF THE METHODS DESCRIBED IN 
FORSYTHE AND MOLFP, "COMPUTER SOLUTION OF LINEAR 
ALGEBRAIC SYSTEMS", PRENTICE HALL, 1967. 

C 

C INDEX OF DOUBLE PRECISION LINEAR SYST.EMS PACKAGE: 

INTEGER MOOE , N, I D IM, I PS < N ) 

RE AL *3 A( IDI M,N), LU U DIM, N I * BIN), XI M I 
REAL DIGITS 

C USES SUBROUTINES DEC MP2 , SOLVE?, AND T MPRV 2 (IE REQUESTED) TO 

C FIND THE SOLUTION TO THE LINFAR SYSTEMS A*X=B. THE MATPIX LU 

C l A NO IPS) WILL CONTAIN THE DECOMPOSITION OF A, 

C THE SOLUTION IS TO DOUBLE PRECISION ACCURACY. 

C MODE = l: DECOMPOSE A, DO NOT IMPROVE THE SOLUTION X. 

C MODE = 2: ASSUME LU CONTAINS THE DEG OMPOSI T I ON, 

C DD NOT IMPROVE X. 

C MODE = 3: DECDMPOSF A AND IMPROVE THE SOLUTION X. 

C MODE = 4: ASSUME LU CONTAINS THE DEG CM POSI T I ON AND IMPROVE X. 

EXTERNAL DEC MP? , SOLVE2, I MPRV2 


SUBROUTINE OEC.M P2 IN, A, I DIM, LU * IPS, *,*) 


INTEGER N, IDIM, IPS(N) 

REALMS At 1 0 1 M , N ) * LU ( I DI M, N ) , DMAX1 , DABS 
C DECOMPOSE THE N*N MATRIX A T MTO TRIANGULAR L £ U SO THAT 

C L*U = A. IPS IS THE POW PIVOT VECTOR. 

C MATPIX A WILL BF OVERWRITTEN BY LU IF A AND LU ARE 

C DECLARED TO BE THE SAME' MAT* IX IN THE CALL OF DECMP2. 

C RETURN l FOR ALL ZERO ELEMF NTS IN A ROW. 

C RETURN 2 FOR ZERO PIVOT. 

C 1 


C 

SUBROUTINE SOLVE? IN, LU, IDIM, B, X, IPS) 


INTEGER N, I D I M , IPS(N) 

RFAL*B LUflOIM.N), BIN), XI N ) 

C SDL V FS A*X*V USING LU FROM SUBROUTINE OECMP2. 

C FI PST SOLVES THE TRIANGULAR LINEAR SYSTEM LY = B 

C AND THFN SOLVES THE SYSTEM U*X = Y. 

C IPS IS THE ROW INTERCHANGE VECTOR FROM DEC MP2. 

1 


C REFERENCES: 
C 

c 

c 

c 


- 200 - 



c +*******«**«************ + ****'** START or PfUTi IPTOTI. *** * ***** * «***' «+ 
f IP TO ?/'<>/ 73 

SUP "OUT 1ST op? ijT ( A , R ) 


C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

G 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c. 

c 

c 

c 

c 


OP PUT# I pt OT L AMI THE FOLLOWING BLOCK O' T A SUBPROGRAM 
A»r USED T 1 GE THF R PQR A CCUN'UL AT I MG EUCLIDEAN INNER PRODUCTS 
OP OOUBLF PRrCl SI ON VECTORS. TMcy MAY ALSO BE USED FOR 
CALCULATING OTHPR EXTPNOFO SUMMATIONS. 


PR OGRANMF R: GORDON GUI LA HORN 

STANFORD UNI VFRSITY 
AUGUST 1G6P 

REVISED BY: mjcuaFL MALCOLM 

COMPUTER SCI EMC F DFPAPTMENT 
STANFORD UNI VPR SITY 
16 NOVEMBER 1^6^ 

p ppgR cf.jcF j t. A nj A. GDR I TH M FOR FLOATING POI A‘T ACCUMULATION 
OF S J MS WITH SMALL RELATIVE ERROR" 

COMMUNICATIONS OF THE ACM, VOL l A NO 11, NOV. 1971 

THE SUM is ACCUMULATED in AN A' R A Y 0 , WHICH IS IN THF 
LABELLED COMMON 1 OP A CCC • • 

THE SUBROUTINES DPERATE AS FOLLOWS: 

BLOCK DATA : INITIALIZES THF ACCUMULATOR *P* TO ZERO. 

THE BLOCK DATA SUBPROGR AM OPE» ATFS AS A 
DATA STATEMENT. THE USER SHOULD NOT ATTEMPT 
TO CALL IT 

DP PUT { A, 0) : ADDS THE PRODUCT A*B T} THE ACCUMUL ATOR 

I P TOTL (X I : COMPUTES THE TOTAL OF THE PRODUCTS IN THF 

ACCUMULATOR AND ASSIGNS THIS VALUE TO X. 

IPTOTI THEN RESETS THE ACCUMULATORS TO 7ERC. 

TP NO OVERFLOWS 1 R UNDERFLOWS OCCUR, THE F ! N A I RESULT IS GUAR- 
A^tfpd T° HAVE 12 S I G'l I P I C AN T DIGITS. IF E! T UE r A DP B IS SMALLER 
» EXPONENT UNDCR PLOW IS LIKELY TO OCCUR IN 
Sfp ABATING THE HIGH AND LOW C B OFR PARTS OF A OR 3.'-. 

SIMILARLY, UNOFPFLOW may OCCUR IF j A*Bj < 1 D**- SO. 


DRKJEJAB PAGE IS 
US POOR jQIJALTTXi 
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C 

C 


SUB POUT I'|f G AL4NC { MM, N, A, LOW, I GH, SCALE) 


INTEGER l , J ,K,L ,M, M, JJ, MM, IGH,LOW« IE.XC 
REAL*8 AfNM.N) , SCALE(N| 

RE4L+P C * F » O , R » S,R2,RADIX 
RFAL* p - PARS 
LOGICAL NOCONV 

THIS SUBROJT I Nr IS A TP ANS L AT ION OF THE ALGOL PROCEDURE BALAMC E * 
NUM. MATH, 11, 293-30'.-! 1969) BY P AR L FT T AMD RFINSCH. 

HANDBOOK TOR AJTO, COMP., V CL . I I -L IN E A* ALGEBRA, 3 1 5-32 6 < 1 971 ) • 

THIS SUBROUTINE BALANCES A REAL MATRIX AND ISOLATFS 
EIGENVALUES WHENEVER POSSIBLE. 

ON INPUT: 

NM MUST BE SET TO the ROW DTMFNSION OF TWO-DIMENSIONAL 
ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM 
DIMENSION STATEMENT; 

N IS THE ORDER OF THE MATRIX; 

A CONTAINS THE INPUT MATRIX TO BE BALANCED. 

ON OUTPUT: 

A CONTAINS T H c BALANCED MATRIX; 

LOW AN D I GH ARE TWO INTEGERS SUCH THAT A(T,J) 

IS EQUAL TO ZERO IF 

(1) I IS GREATER THAN J AND 

(2) J=l,. LOW-1 OR I=IGHH,..,»N; 

SCALE CONTAINS INFORMATION DETERMINING THE 
PERMUTATIONS AND SCALING FACTORS USED. 

SUPPOSE THAT THE PRINCIPAL SUBMATPIX IN ROWS LOW THROUGH I GH 
HAS BEEN BALANCED, THAT P ! J ) DFMOTES HE IND C X INTERCHANGED 
WITH J DURING THE PERMUTATION STEP, AND THAT THE FLF^MTS 
OF TH f DIAGONAL MATRIX USED ARE DENOTED BY D( I , J). THEN 
SCALE! J) - D ( J ) , FOR J = 1 , • • • , LD W- l 
D ( J, J) , . . ; LOW,..., I GH 

- P ( J ) J = IGH+1 , N. 

THE ORDER IN WHICH THE INTERCHANGES ARE MADE IS N TO IGH+1, 

THEN 1 TO LOW-1 . 

NOTE THAT 1 IS RETURNED FC» IGH IF I GH IS ZERO FORMALLY. 

THE ALG^L PROCEDURE EXC CONTAINED IN BALANCE APPEARS IN 
BALANC IN L INE * (NOTE THA T THE ALGOL POLES OF I D F N T I F I F P S 
K»L HAVF BEEN R C VERS C D. ) 

QUESTIONS AND COMMENTS SHOUL 0 BE DIRECTED TO r, c. GA’BDWt 
APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LAB OP A TORY 
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c 

c 



SUBROUTINE ELMHES ( NM,N, LOW, IGH» At I NT) 


INTEGER I , J . M,N , LA ,NM, I GH, KP1 , LOW ,MM1 ,MPl 

REAL *8 A(NM.N) 

REALMS X » Y 
RFAL*8 DABS 
INTFGER INT(IGH) 


THIS SUBROUTINE IS A TRANSLATION OF THE 
NUM. MATH. 12, 349-368(19681 BY MARTIN 
HANDBOOK FOP AJTO. COMP., V OL. I I-L INE AR 


ALGOL PROCEDURE ELMHES 
AND WILKINSON. 

ALGEBRA, 339-3581 1971 ) 



GIVEN A REAL GENERAL MATRIX, THIS SUBRDUT INE 
r cqu:FS A 5U8 M A to I X SITUATED IN ROWS AND COLUMNS 
LOW THROUGH IGH TO UPPER HESSENBERG FORM BY 
STABILIZED ELEMENTARY SIMILARITY TP ANSF OPM ATI QMS • 


ON INPUT: 

NM MUST BE SET TO THE ROW DIMENSION OF TW3-D I MFNSI DNAL 
ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM 
DIMENSION STATEMENT; 

N IS THE ORDER OF THE MATRIX; 

LOW ANO IGH ARE INTEGERS DETERMINED BY THE BALANCING 
SUBROUTINE BALANC. IF 8 AL ANC HAS NOT 8EFN USEO, 
SET LOW=l, igh=n; 

A CONTAINS T H p INPUT MATRIX. 


ON OUTPUT: 

A COMfAINS THE HE S SENPER G M ATR IX . T HF MULTIPLIERS 
WHICH W E° E USED IN the REDUCTION ARE STORED IN THE 
REMAINING TRIANGLE UNDER THE HESSENBERG MATRIX; 

INT CONTAINS INFORMATION ON THE ROWS AND COLUMNS 
IN T EPCHAN3ED IN THE REDUCTION. 

ONLY ELEMENTS LOW THROUGH IGH ARE USED. 

quc$t ions and comments should be directed to b. s. gar bow, 

APPLIED MATHEMATICS DIVISION, ARGQNNE NATIONAL LABORATORY 
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r. 

c 

c 

c 

c _ 

c 

c 

C~ 

c 

c 

c 

c 

c 

c 

c 

c 

G 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


SUB D DUT 1 ME HQ3 { N v , N, L CW 1 1 GH ,H • WR , W 1 v IERP ) 

I NT ^or R ! , J ,K,Lt M, N, EN* LL » MM. N A .NM , IGH, ITS. tnw,MP2» EM*? » I ERR 
RFAl*° H( MM ,m )» WR I N5 ,WI IN) 

RFAL*R P, Q » R »S »T , W .X.Y. ZZ» MAC HEP 
REM*P DSCRT.DABS.QSIGN 
INTEGER NINO 
LOGICAL NQTLAS 

THIS SUBROUTINE IS A TRANSLATION OF THF ALGOL PROCEDURE HQR, 
NU W . m ATh. 1A, 2 1 9-2 3 1 ( 1970 ) BY MARTIN, PETE C S, AND WILKINSON. 
HANDBOOK POP AJTO. COMP,. VOL • T I-L INEA * ALGEBRA, 359-3711 19711. 

THIS SUBROUTINE FINDS THE E IGFNV , U^S DF A RpAL 
UPPER HESSE NB ERG MATRIX BY THE QR METHOD* 

ON INPUTi 

NM M«JST 9T SET TO THE ROW DIMENSION OF TWD-OI MPNSI DNAL 
ARRAY PARAMETERS AS DECLARED IN THF CALLING PROGRAM 
DIMENSION STATEMENT*. 

N IS THE ORDER O e THE MATRIX; 

LOW AND IGH AO F INTEGERS DETERMINED BY THE BALANCING 
SUBROUTINE BALANC. IF B AL ANC HAS NOT BEEN USED. 

SET L0W=1, IGH=N; 

H CONTAINS THE UPPER HESSENBERG MATRIX. INFORMATION ABOUT 
THE TRANSFORMATIONS USED IN THE REDUCTION TO HESSENBERG 
FORM BY ELMHFS OR ORTHFS, IF PERFORMED, IS STORED 
IN THE REMAINING TRIANGLE UNDER THE HESS ENBFPG MATRIX. 

ON OUTPUT : 

H HAS BEEN DESTROYED. THEREFORE, IT MUST BE SAVED 
BEFORE CALLING HQR IF SUBSEQUENT CALCULATION AND 
BACK TRANSFORMATION OF EIGENVECTORS IS TO BF PERFORMED; 

WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS., 

RESPECTIVELY, OF THE EIGENVALUES. THF EIGFMVALUES 
APE UNORDERED E XCF PT THAT COMPLEX CONJUGATE PAIRS 
OF VALUES APPEAR CONSECUTIVELY WITH THE EIGENVALUE 
HAVING THE POSITIVE IMAGINARY PART FIRST, IF AN 
FRODR EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT 
EC. R INDICES IERR+l .... ,N; 

IF PR I S SET TO 

ZERO FOR NORMAL RETURN, 

j IF THE J-TH EIGENVALUE HAS NOT BEEN 

DETERMINED AFTER 30 ITERATIONS. 

QUESTIONS AND COMMENTS SHOULD B F DIRECTED TO B • $• GAPBOW* 
APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY 
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o p n o o o o n o o o o o o o o o o o o o o o n o o n r> o o ooo 


OP IV PIG RC'JTINE EPR UNCONSTRA TNEO MINIMIZATION U ! TH- n UT DERIVATIVES 


SURRO l J T INP MINM7.FI N » NLDIM, FUN, S'TART ) 


IMPLICIT PFAM'8 ( A-H, L, O-l ) 

LOGICAL P 0 IFF , tr-NV 

C CM MOM /GFITFP/ ^TJ, TOL, STFPMX 
COMMON /LOG! CL/ BRIEF 

DIMENSION X I ?0 ) , D ( ?0 ) • LI450), HI30I 

EXTERNAL FUN, START 


THIS SUBROUTINE USES SUBROUTINE CNMOIF TO MINIMIZE 4 FUNCTION 
OF N VARIABLES, WITHOUT USING DFFIVATIVFS. 

THE CALLING PROGRAM MUST DECLARE FXTERNAL THE SUBROUTINES EQUI- 
VALENT TO FUN AND START. THE CALL IS THEN OF THF FOPM 
CALL MINMZEI N, NLDIM, FUN, START ). 

PROGRAMMER: DAVID SAUNDERS, STANFORD UNIVERSITY. 

LAST UPDATED: APRIL, \973. 

OTHER ROUTINES REQUIRED: 

, QNMOI'F, LINSCH, MCDCHL , APROXG, 

OUTPUT, FUN, START. 

STORAGE FOR ARRAYS X, D I DIMENSION NI AND L (DIMENSION NLDl M = 
M# ( N- I ) /? I may RE SUPPLIED BY THE CALLER. HOWFVER, IT WAS CON- 
SIDERED CLEANER TO DIMENSION ABSPLUTELY ANY LOCAL ARRAYS IN THE 
OTHER SUBROUTM ES , NO MULT I-D I ME N SI ON AL ARRAYS APE USED SO THE 
USER HAS ONLY TO GIVE ALL LOCAL ARRAYS ABSOLUTE DIMENSIONS OF 
(AT LFAST) N, WITH ONE EXCEPTION: 

MQDCHL USES ARRAY S WITH DIMENSION N+l • 

THE USER MUST PROVIDE HIS OWN EQUIVALENTS OF THE 1 SU ac OUTINFS 
FUN AND START SUPPLIED HERE AS ILLUSTRATION. 

COMPILERS: WAT FI V OP 0S/360 FORTRAN H I PPT=2 0 EC n MMENPF D I 



of*oooo'of'ooor> 


C 


c- 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c' 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


SUBROUT I ME ON MO I F ( N» NLDIM, F, X, L, 0, H, CONV, FUM ) 


IMPLICIT REALMS (A-H, L» O-l > 

logical umtl, print, pptff, conv, succes 

COMMON /CRITFP/ FT A , TOL, STEPMX, OEPS 

COMMON A'UM8«$/ NtJMF 

COMMON /COUNTS/ NFEVAL, ICOUNT 

COMMON /TOTALS/ MFTOTL, NITERS 

COMMON /LOGICL/ 8P l F F , PRINT, UNI TE 

0 I MPN SI ON X(N), LINLOIMI, DIN), -UN) 

DIMENSION GK (30 ) , GKPLS1I30I, W(30), POO), PPOO) 

EXTERNAL FUN 


THIS SUBROUTINE ATTEMPTS UNCONSTRAINED FUNCTION MINIMIZATION, 
USING A REVISED QI.IAS I -NEWTON METHOD WITHOUT DERIVATIVES. 

IT ISA TRANSLATION QF THF ALGOL PROCEDURE ON MO' IFF IN THE RE- 
PORT BELOW, TO WHICH THE USER IS REFERRED: 

R frp£ R eNCF : "IMPLEMENTATION OF TWO REVISED QUASI-NEWTON 

ALGORITHMS FOR UNCONSTRAINED OPTIMIZATION." 

BY GILL, MURRAY, AND P ITFIFLD (APRIL 1972). 
(REPORT DNAC II OF THE NPL, LONDON.) 

THE ROUTINE SEEKS THE POINT X AT WHICH THE TWICE CONTINUOUSLY 
DIFFER ENT I ABLE FUNCTION F ( X ) ATTAINS ITS LEAST VALUF. IDEALLY 
THF VARIABLES SHOULD BE SCALFn SO THAT THE HESSIAN MATRIX OF 
SECOND DERIVATIVES AT THE SOLUTION IS APPROXIMATELY ROW-EQUI- 
LIBRATED, WITH THE FUNCTION MULTIPLIED PY A SCALAR SUCH THAT 
IT ACHIEVES A MAXI M U M VALUE OF UNITY WITHIN A UNIT SPHERE A- 
ROUNO THE MINIMUM. IT MAY NOT BE POSSIBLE TO FULFILL EITHFP^ 
OF THFS c REQUIREMENTS. 

GIV^M AN INITIAL APPROXIMATION TO THE MINIMUM AND AN ESTIMATE 
( LOW c R BOUND) 3E THE MINIMUM VALUE, THE ROUTINE CALCULATES A 
LOWER FUNCTION VALUE AT EACH ITERATION. WHEN THF CONV pp GFNCF 
CRITERIA APE SATISFIED. THE ROUTINE GIVES THE ESTIMATED POSI- 
TION OF THF MINIMUM, THE FINAL FUNCTION VALUE, AND THE FINAL 
CHOLES KY F ACTOR I ZATIDN OF THE APPROXIMATE HESSIAN MATRIX, 

rw c p £p »mc J r.? 5 AND CDMMnM VARIABLES OF SUBROUTINE QNMQIF APE 
INI T I A L I ZED IN SUBROUTINE STAPT. A D " SC R I p T I ON OF THESE IS 
GIVEN NOW: 




INPUT TO SUBROUTINE QNMDIF: 


K| : the \‘!J M B F R OF V AR I A PL F S X ( I ) OF THE P UNCTIDN P(X). 

pyvj . a SUB = 0 U T I N c RETURNING T HE VALUE DP THE C UNC T I ON , AT 

THE POINT D c F!NfD BY VECTOR X(I). THE CALL IS OF THE 
F T P v CAL L FUN ( N, X , F ) . 

p • A L n WE R GO UN D FOR THE MIN I MUM V *1. UE OF THE FUNCTI ON* 

X : AN INITIAL ESTIMATE OF THE SOLUTION, 



OOOUOOOUOOUUOUOUUOOOOOOOOO o<-><-HJuOOOOUUOUOOUOOUOOOU 


L. D : ACRAYS WHICH WILL NORMALLY NOT BE INITIALIZED BY THP 

t j s: r ^ i-vut s fp "unitl" r^low). 

H • : A 1 <>'‘1 4 '■) 3 4 Y CONTAINING THF INTERVALS FOR P I FFF RF NC l b G 

F(X) ALONG PACH OF THF COORD INATE DIRECTIONS. TYPICAL 
V Al (IF FOQ THE HI!) IS DSQPTI DE PS )• 

DEPS : THE D EL AT I VE MACHINE P BC C ISION. 

pTA : THE T p L M I A'AT I Of.' C P I TE° I ON FOR T H c LINFAP SEARCH. IT 

SHOIIl D MAV C A VALUE IN THF RANCF 0 TO 1. THE CLOSER 

TO Z cp D IT IS. THE GRF ATER THE NUMBER OF EVALUATIONS 

OF THE Fit NOT ION THAT WILL BE PERFORMED, WHILE THE 
CLOSER Tn i [T IS. THF GREATER THE NUMBER OF ITERA- 
TIONS LIKELY. ETA = 0.1 IS SUGGESTED TO START WITH. 
TOL : T HF OVERALL TERMINATION CRITERION (NORM OF THF GRADI- 
ENT). A TYPICAL VALUE IS 1C-6. (A GOOD FSTIMATR IS 
APPROXIMATELY DSC 0 T ( DEPS ), BJT THIS CAN BE p Fl AXED 
(INCRFASRD) IF FEWEP SIGNIFICANT FIGURES ARE ACCEPT- 

A BL F . ) . 

STFPNXS AN UPP=R BOUND ON TH p STFP ALLOWED ALONG A DI oc C.TICN 
IF S p A*C-H. THIS CAM BE USED TO PPEV C NT OVERFLOW IN 
IN THE COMPUTATION. IF AN APPROXIMATE SOLUTION ISN'T 
KNOWN, AND OVERFLOW IN COMPUTING THE FUNCTION IS UN- 
t IK p LY, THEN RTFPMX CAN BE SET VERY LARGE (!Dll. SAY) 
SO THAT IT WILL NOT INFLUENCE THE ALGORITHM AT ALL. 
IF THE SOLUTION IS KNOWN TO BE WITHIN A CFRTAIN RANGE 
OF THE INITIAL ESTIVATE, THEN STFPMX CAN BE SUITABLY 
LOWERED • 


C 


UNI TL : A LOGICAL VARIABLE WHICH SHOULD BF SET TO TRUE UNLESS I 

AM APPROXIMATION OTHE c THAN THE UNIT MATRIX IS KNOWN I 

FOR the HESSIAN MATRIX OF SECOND DERIVATIVES OF FIX). I 
’IF UNITE IS TRUE, THE UNIT MATRIX IS SUPPLIED BY THE ) 
SUBROUTINE, ELSE IT IS ASSUMED THAT AN LDt ' FACTOPIZ- I 
AT ION IS G I V r M IN TH° ARRAYS L AND 0. I 

PRINT : S FT THIS L n GI CAL VARIABLE TO TRUE IF COMPLETE OUTPUT I 

IS POSITED APTER FVERY ITERATION. S 

BRIEF : S FT ^HIS TQ TOIF IF MPR c CONCISF OUTPUT IS PR FFF p RFC. | 

FULL INFORMATION FROM T HE FINAL I T c RA T ION ' IS STILL I 

ooo VIDEO EVEN IF PRINT IS FALSE. I 


OUTPUT FROM SUB P DU TINE GNMDIF: 

NFTOTl : THF TOTAL NJMBFR OF FUNCTION EVALUATIONS USED, 

NITERS: THE TOT AL NU‘-‘REP OF L I NEAR SEARCHES PERFORMED. 

C n NV : A LOGICAL VARIABLE SET TO TRUE IF TF P fl I NATION OCCUFS 
WITH THF CONVERGENCE CRITERIA SATISFIED, AND EAI.SF IF 
A LOWER POINT CANNOT BE FOUND ALONG A PARTICULAR OIP- 
f r T ION OF SFA°C H. 

ALSO OUTPUT ARE FINAL VALUES F0° F. X, P. AND L. 
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APPENDIX C 


HYDRAULIC TORQUER DYNAMICS 

A schematic of the torquer is shown in Fig. C-l. It consists of 
two parts: 


(a) A hydraulic rotary actuator mounted on the gimbal, 


(b) A twor-stage electrohydraulic valve that controls the fluid 
flow to the actuator. Electrohydraulic valves of this type 
are described in the literature (e.g., GUI-l). The output 
torque of the actuator is proportional to the differential 
pressure across its vanes.' 

The first stage of the electrohydraulic valve is a nozzle and flapper 
driven by a d-c torque motor. Its output differential pressure, 

Pcq " p c 2 > "‘• s proportional to the torque motor current. 

The second stage is a two-orifice spool that is displaced by the 
applied and feedback differential pressures. The fluid flow to the 
actuator is proportional to the displacement. 

The dynamic equations of the torquer are given below. For the 
spool. 


m x - p A, - pA - bx , 
v v ■ c 1 2 v ’ 


(C-l) 


where 


m = 
v 

P c = 
P = 


A. , A 0 
lv 2 


b = 


spool displacement 
spool mass : 

P'd _ p C2i t P e command differential pressure 
p^ - P^, the feedback differential pressure 

the areas of application of the command and feedback pressures 
viscous friction coefficient. 
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Assuming a first order lag between the flow rate and the displace- 
ment, we get 

4 = c x - — q , (C-2) 

T 

V 

where 

q = load flow rate 

c.^ = proportionality constant. 

Neglecting the mass of the spool in Eq. (C-l) and assuming A = A =; A 

X u 

(required for p = p^ in the steady state), the final equation for the 
spool is 

<1 = K A ( P C - p) - ± q ( 0 - 3 ) 

V 


where 



For the actuator, 

P = ~ (q - - V), '■= (C-4) 

e 

where 

B = bulk modules of the fluid 

v = volume of the actuator chamber 

. e ■■■■■• 

k^ = leakage coefficient across the vanes 
D = chamber volume change per unit angular displacement 

Pi 

6 = relative angular displacement between the stationary and 

moving parts of the actuator. 

Equations ( C-3) and (C-4) are the state equations for the torquer 
with p and q as states* They are rewritten below 


- 210 - 



- B -*l 
v L 

e 


- k. 


B_ 

v 


v-* 


0 


P c + 


B D. 


v 


0 0 


The transfer functions from and 0 to the output p: 


P_ 

P, 


k A (B/v e ) 


^ ^ (f x„ * f)‘ ' f ^ * *a) 


N x (s) 

dTsT 


(c-5) 


Numerical values: 


P 

8 


V T 
e v 


V s) 


d(s) D(s) 


(C-6) 


B_ 

v 


= 2.75 X 10 


5 lb 


xn 


k_ = 5 X 10 


. . 5 

-4 in 


lb sec 


\ = 1 - 4 


. 5 
in 


lb sec 


1 -1 
— = 24 sec 


■ 3 / - 

D = 0.3 in /rad 
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1 



The damping of this second order system depends mainly on the leakage 
across the vanes and therefore cannot be determined precisely. The 
natural frequency is proportional to the bulk modulus which may vary 
considerably with environmental conditions. 

To use the torquer dynamic equations in the system, it is convenient 
to replace the command and output pressures by command and output accel- 
erations. The relationship is 



(C-7) 


Where 

a = acceleration 
I = load moment of inertia. 

Substituting (C-7) into (C-3) and (C-4), we get the state equations 



For the inner azimuth gimbal ( I = 


in-lb-sec 
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. . 2v 

For the elevation gimbal ( I = I g = 20 in -lb-sec 



where cp and £ are defined in Figure V-l. 




APPENDIX D 


NUMERICAL DATA FOR THE TRACKING TELESCOPE 


ik 

1 . Moments of Inertia 

Telescope plus inner azimuth gimbal about lb axis: 

2 

I = 680 in-lb-sec . 

Telescope plus inner azimuth and elevation gimbals about 2b 
axis: 

2 

I = 400 in-lb-sec . 

2 


Telescope plus inner azimuth and elevation gimbals about 3b 
or lb axes: 

2 

I = 920 in- lb-sec 

Outer azimuth gimbal about 1 axis: 

2 

I = 2500 in-lb-sec 

4 

2. Acceleration Limits 

,2 

Inner azimuth: 4 rad/sec 

> 2 

Elevation: 3.5 rad/ sec 

/ 2 

Outer azimuth: 0.65 rad/sec . 


3 


Measurements 


(a) Target detector 


Noise: 100 (irad rms at 2000 sec 

in -U ,2 

r, = 10 rad sec 
d 



See Fig. V-l for definition of axes. 
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(b) Gyro 

Time constant: 1,5 msec 

Gain : 2.5 

Noise: 5 prad rms at 3000 sec 

r : 2 X 10 -1 '* rad 2 sec 

g 

4 . Disturbances 

Disturbance torques on inner gimbal: 


2 —I 

a T = 0.1 rad/sec at 200 sec 


-4 2 3 

10 rad /sec 


Torquer noise (1% of full scale output) 


/ 2 -1 
0.04 rad/sec at 300 sec 


N 

, n -5 .2/ 3 

q = 10 rad /sec 


5. Angular Limits 


Cp = + 3 C 

T max — 


g = -6 Q , + 25' 

max 


unlimited. 
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APPENDIX E 


STRAIGHT LINE FLYBY 


The geometry of a straight line flyby is shown in Fig. E-l below 



FIG. E-l GEOMETRY OF A STRAIGHT LINE FLYBY 


We have 


tan 0 



where t is measured from the time of passage at 0. 


e = 


V cos 0 
P 

„2 


V 2 „ VI 

oos 0 = R -— 2 

1 + T 


0 = 2 — jt cos 0 sin 0 = 2 — rj 

R 2 R 2 


V 2 T 


,, 2 2 
(1 + T ) 


' v 1 

0 = 0. 65 — - at t = — • 

max R 2 JT 


0 as a function of t is shown in Fig. E-2. 


The maximum angular acceleration in aximuth that the system can 
accommodate before the torquer of the outer azimuth gimbal saturates 
is 0.65 rad/sec . For this acceleration we obtain 
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The time at which the maximum acceleration occurs for this ratio is 



0.58. 


The rate of change of the acceleration is low relative to the expected 
system bandwidth. If a Type 2 (constant acceleration error) system 
is specified, the error will therefore be approximately proportional 
to the instantaneous acceleration. The acceptable error at the point 
of maximum acceleration (t = l//3) can therefore be specified as the 
maximum error for a constant angular acceleration command. This 
error should be of the order of the measurement noise. 




APPENDIX F 


SYSTEM DYNAMIC EQUATIONS 


The dynamic equations of the system are found in this Appendix 
using the Lagrangian method. Referring to Fig. V-l , and to the moments 
of inertia defined in Appendix D ; the kinetic energy of the system is 


T = I 1 ( cp + -f cos e ) 2 + ^(l 3 - cos 2 e + 


(F-l) 


1*2 1-2 .2 1 *2 
+ - I 3 e + - * sin 6- + 2 y • 


The potential energy is V = 0; .*. L = T. 

The moments about the gimbal axes are given by 




M 

e 



Sr 

Sep 

Sr 

"Sep 

St 

Se 


_d_ St St 

dt Sep Sep 

d_ St St 

dt Sg Sg 

d_ St _ St 
dt Sjf cty 

• • " 

I (cp + \|r cos e) 
0 


- 219 - 


£ &\st 


uT *2 • • • 

— = (I - I + I )\jr cos e sin e - I (cp + \|r cos s)\jr sin e 

^ o X 1 

oe 


• • • 2 * 2 * 

= I (cp + ijr cos e)cos e + (I - I, )i}r cos e + I_ \Jr sin g + I Jf 
1 3 1 2 4 


= 0 . 


Defining a, = cp + ijr cos e as the total rate about the inner 
azimuth axis, the moment equations are: 


M 


9 


la 


M 


M 


I e ■+ I^cty sin e -(I g ^ - I 3 H cos e sin e 

2 2 •• 

[ (I 3 - I 1 )cos e + I 2 Sin e + I 4-Jt + I 1 a cos e 
+ 2(1^ - I3 + 6 cos € sin e - I^Gle sin e . 



APPENDIX G 


LINEARIZATION OE THE ZERO ORDER HOLD (ZOH) 


The transfer function of the ZOH (zero order hold) is [OG-l] 


G (s) = 

h s 


-Ts 


(G-l) 


The Laplace transform of the sampled signal x (t) is [OG-l] 

00 

(s) = i ^ x(s + k) . (G-2) 


x* 1 


k=-“ 


If most of the incoming signal energy is in frequencies below the 
sampling frequency, x*(s) can be approximated by 


x*(s) « - x(s). 


(G-3) 


Combining Eqs. (G-l) and (G-3), the transfer function for the sampler 
and hold is 


G (s) = 

sh 


1 - e 


-Ts 


Ts 


1 - e 


a 


-a 


where 


a = Ts = 


2:rts 


(G-4) 


Using 


. , a I, a/2 -a/ 2 S 

sin h - = g (e - e ) , 


(G-4) may be rewritten as 


sin h 


G sh <s) 


a/2 


a 

2 -a /2 

—■ e 


(G-5 ) 
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Two approximations to 


G , (s) will now be examined 
sh 


-a 


(a) Approximation of e in Eq. (G-4) by 


-a „ i - a/2 

e — 5 


i + a / 2 


—a/ 2 

(b) Approximation of e ' in Eq. (G-5) by 

-a/ 2 l - a/4 

G — 1 • 

1 + a/4 


The transfer functions for these approximations are 


G (s) = 

a 


TtS 

1 + to 


0 


G b ( S ) = 


, its , its 
sxn h — 1 - - — 

co 2co 

0 0 


its 

to„ 


1 + 


TtS 

2to, 


(G-6) 


(G-7) 


For s = jco, the amplitudes and phases for these approximations are 
given in Table G-l. 


Table G-l 



Exact 

Approx, a 

Approx, b 

A 

. TfCO 

sm — 
CO 

0 

• 1 

| 

Exact 

TtCO 

CO 

0 


> 

' Cx) 

Jt 

“o 

■ -1 TtCO 

tan — 

,J) o 

2 t a „- 1 = - 

: 2 “ 0 : ;; 


- 222 - 




Numerical values for the amplitudes and phases at different fre- 
quencies are given in Table G-2. 


Table G-2 



4 l 

91 ) 

CO 

Exact 

Approx. 

Approx. 


Approx. 

Approx. 

w 0 


a 

b 

Exact 

a 

b 

0.125 

0.974 

0.93 

0.974 

22.5 

mm 

22.2 

0.25 

0.9 

0.787 

0.9 

45.0 


42.8 

0.5 

0.63 

0.53 

0.63 

90.0 

m 

76.0 


From Table G-2 it can be seen that if approximation b is used with 
a fixed gain of 0.9, the amplitude and phase correspond quite closely 
to the exact system up to w = 0.25 co . The transfer function for 
this approximation can be written as 


G b (s> 


0.9 


4£ o 

4f o 


- S 


+ S 


(G-8) 


where f^ is the sampling frequency. 

The block diagram of this transfer function is shown in Fig. G-l, 
Its state representation is obtained from the block diagram as 


P = - 4f 6 -7.2f a + 7 . 2 fa 

0 0 0 c 

y = p + 0.9a - 0.9a c 


(G-9) 
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APPENDIX H 


GYRO DYNAMICS 


The dynamic equation of a rate integrating gyro (RIG), which 
represents the sum of the moments about its output axis, is 


J0 O = - O0 o + Hu t + 


where 


0 

c = 


0_ = output axis deflection 

damping coefficient about the output axis 

H = angular momentum of the gyro 

to = angular rate about the input axis 

k^ = torquer scale factor 

i^ = torquer current 

J = inertia about the output axis. 

For a typical RIG (e.g., Honeywell GG 334), the numerical values of 
the parameters are 

4 2 

H = 5 X 10 gm-cm / sec 

4 2 , 

c = 2X10 gm-cm / sec 

2 

J = 30 gm-cm 

2/2 

= 237 gm-cm /(sec ma) . 

The transfer function of the gyro is therefore 


s(s + j) 


H T . 

, — W. + — - i 

c. VJ l J 


- 225 - 



With the given values of the parameters 


e - [1660 CO. + 7.9 i]. 

0 s (s + 660) i 

If c/j is much larger than the system bandwidth, the transfer function 
from the input axis rate to the output axis deflection is approximately 


V s) 


H 


( 0 . (s) Js 
1 


If the torquer current is made proportional to the output angle, 


1 = ~ Hk A 9 0 ’ 


the gyro becomes a rate gyro. For this case, 


JS 0 = -ce o . „k A k T e 0 + Ho.. , 


and the transfer function is 


V s) 

CO. (s) 

1 


H 1 

J 2 c H , ' 

S + J S + J k A k T 


The output angle is considered as a measure of the input axis rate, 
The steady state scale factor is 


e 


Oss 


CO 


xss 


k k * 
A T 


In practice, the torquer current and not the output axis angle is meas- 
ured. . 

In control applications the dynamics of the rate gyro is often 
neglected and the output angle is considered as a direct measurement 
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of the input rate, viz., 


V s) 


k A k T 


w.(s) . 

X 


The influence of this neglection is to couple the system and gyro 
eigenvalues with the coupling becoming less pronounced as the gain 
of the rate gyro is increased. 

For the system using the reduced order estimator, a root locus 
as a function of the gyro gain is shown in Fig. H-l for values of this gain 

’ g g _ 2 

[ (H/j)k Lj ranging from 10 to 10 sec , which corresponds to gyro 
** 3 4 “1 , 7 -2 

natural frequencies of 10 to 10 sec . For gains above 10 sec 

the root movement is hardly perceptible and even for [(H/jJk^k,^] = 10 , 

it is not excessive. 


The influence of this gain on the rate gyro output noise will now 
be considered. The gyro pickoff noise is considered as having an rms of 
5 [xrad at 3000 sec 1 (see Appendix D). The gyro with its pickoff 
noise may therefore be represented by the block diagram in Fig. H-2. 
From this block diagram, the state equation of the gyro and its pickoff 
noise can be represented in state form as 




u> 0 = [k, k, 0]x, 
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fig. h-i root locus as a function of the rate gyro gain 





FIG. H-2 BLOCK DIAGRAM OF RATE GYRO WITH NOISE 














APPENDIX 


THE ROOT SQUARE LOCUS METHOD FOR DETERMINING STATE WEIGHTS 


The closed loop eigenvalues for a system designed by quadratic 
synthesis are the left half plane eigenvalues of (Eq. 2,6) 

B + Y T (-s)AY(s) = O, (^-l) 


where 

Y(s) = (si - f) G 

is the transfer matrix from the control to the state: 

x(s) = Y'( s ) u( s) . 

For single input systems, Eq. (c9-l) can be written as 

b + tr[AY T ( - s)y( s) ] . C' 9 " 2 ) 


For a matrix A with diagonal terms only, this equation becomes 

n 

1+ 2 y i (-s)y.(s) , ( c9 “ 3 ) 

x=l b i i 


where 


; n.(s) 
y ( s ) = 

D( s) 


is the transfer function from u to x.. 
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If the ratios a^/b are fixed for all i ^ k, the root locus 
as a function of a kk / b can be found, by writing Eq. (J-3) as 

D(-s) D(s) + 

+ if V- S) V S > = 0 

or, as 

D 1 (-s)D 1 (s) = - ~ Y k (-s)Y k (s) . (c9-4) 

The eigenvalues of the left and the right hand side are symmetric 
about both the real and the imaginary axis. 


Y.(-s)Y i (s) 
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APPENDIX J 

REDUCED ORDER OBSERVER 


The concept of the reduced order observer was introduced by Luen- 
berger [LU-l], According to this concept, the state of a system can be 
observed by an observer of order n - m where n is the dimension of 
the system and m the number of measurements. 

Various approaches are available for the actual design of the 
observer. The approach used here is based on a design method proposed 
by Gopinath [GO-1], 

Consider the observable system: 


x 0 = Vo + V + V 


y = V + v • 


(J-la) 


A state transformation is performed such that the state vector has the 
form 


and 



y 


= X 

n 


+ v. 



(J-lb) 


In this representation, only the partial state vector has to 

be observed since x is measured directly. 

n 

The dynamic equation of the transformed system can be put in the 

form 
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(J-2) 
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n 


from which the equations for x and x are 

r n 


x = Fx + Bx + G u + T w 
r r r r n r r 


x 


Hx + Fx + G u + T w .. 
r r n n n n 


(J-3a) 

(J-3b) 


From Eq. (J-3b), a "measurement" for the state x can be defined as 

r 


y = y-Fy-Gu = Hx + v 
r n n r r r 


(J-4) 


where 


v = r w + v - F v , 
r n n 


(J-5) 


An observer for x using this "measurement" has the form 
r 


Fx + By + gu + K(y - Hx), 
rr r r rr rr 


U-6) 


Substituting for y from Eq. (J-4), the equation of the observer be- 


comes 


* 

x 


= (F -K H ) x + (B -IC F )y + (G -K G )u + K y . (J“7) 

r r r r r r n r r u r 


In practice, it is of course not feasible to differentiate the 
measurement y and therefore a new state variable is defined as 


6 = x - K y . 

r r 


(J-8) 


The dynamic equation for 6 is 


(F - K H )x + (B - K F )y + (G - K G )u 
r r r r r r n r r n 


(J-9) 


or 
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0 













The error state equation can be obtained by substracting Eq. (J-6) 
(J-3a), and using (J-la). It is 


x = (f - K H ) x - Bv-Kv + P w , 
r r rrr r rr r 


Substituting for from Eq. (J-5a), 


x = (F-KH)S -(B -K P ) V' + (r -K T )w - K v” 
r 'rrr' r rr n rrn' r 


In order to avoid the use of v, a new state variable is defined 


8 = x + K v . 

r r r 


Its dynamic equation is 

• 

e = (f - k h )e -[b - kf +(f - k h )k ]v 

r v r rr'r r rn'r rr'r 

+ (r - k r )w . 

v r r n' 

The control equation, (J-lOa), can now be rewritten as 


u = -C x + c 6 - C X -(c K + C )v . 

rr rr yn'rr y 


Substituting Eq. (J-14) in (J-3a) , and adding (J-13), the dynamic 
tions for the augmented system are obtained. They are 


X 


F -G C 
rrr 

I 

1 

1 

-L. 

B -G C 
r r n 

1 

1 

1 

-i.. 

G C 
r r 

A 

n 


H -G C 
rnr 

1 

1 

1 

1 

1 

F iG C 
n n n 

1 

1 

1 

G C 
n r 

e 

r 


o 

T 

i 

i 

0 

“T 

l 

l 

7 -K H 
r r 


r 


L J 


r -k r 

r r n 


-M. 


-M. 




M. 


, = G (C K + C ) 

1 rrr y 


M. = G (C K + C ) 

2 n r r y 

M = B - K F +(F - K H )K . 

3 r rn r rr r 
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f rom 


(J-ll) 


(J-12) 

as 


(J-13) 


(J-14) 

equa- 


(J-15 ) 


where 



To apply this method to the inner gimbal system as represented in Fig. 
V-2a, the detector measurement e must be defined as a state. 

From Appendix G, 


0 = p + 0.9 a 


$ = -4f 0 P - 7.2f Q a . 


(J-17a) 

(J-17b) 


Differentiating (J-I7a), inserting (J-17b), and using a = w, the state 
equation for e is 


€ = -4f e - 3.6f Q 0i + 0.9w . 


(J-18) 


Using 


x r = [a, a, r] 


< = Ce, «] 


the dynamic equation, (J-2), can be written in the explicit form 
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•md 



(J-19) 


From this equation, the augmented state equation can be set up. The 
estimator gain matrix K is then found by parameter optimization, using 
the PAROPT program. The resulting K matrix is 


K = 


-0.139 

1.03X10 4 

4 

2 

1.46X10 

-8.75X10 

7 

3 

-5.05X10 

2.73X10 
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The observer eigenvalues are: - 36 + 610 j , -0.54. Note the very 

low real eigenvalue. 

These results were used in Chapter V-D-4 for comparative evalu 

ation. 
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APPENDIX K 


DETERMINATION OF ACCELERATION ERROR 

In this Appendix, the acceleration errors are determined for two 
types of systems; (1) Type DR I, (2) Type DR H2. The errors are deter- 
mined by computing the value of the integral state i for a constant 
rate command. This is equal to the value of the state e = di/dt for 
a constant acceleration command. For a constant rate command, g = 0 
(Type II system). 

1. Type DR I 

No control is required for constant rate and the sum of the inputs 
to the gyro must therefore be zero. From Fig. V-4 this sum is 

£ u = a i + a_ e + w = 0. 

0 1c 

Since g = 0 and to to 

c 

This is also the value of the error e for a constant acceleration 
command and the acceleration error coefficient is therefore 

e _ 1 

a ~ K, 

c i 

2. Type DR H2 (Fig. V-6c) 

The sum of the controls that are obtained from the nonzero estimator 
states and the integral control must be zero, 

From Fig. V-7b, it is obvious that for constant to 
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/s 

OJ 


0) 


From Fig. V-7a, the state equation for estimator 1 is 




- 



- - 


r ~ 


n 

A 

P 


-4f -k 
0 1 

-7. 2f -0.9k, 
0 1 


A 

P 
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k 

1 

e + 

0 

A 

a_ 


- " k 2 

-°.9k 2 . 


A 

a 




1 


A /\ 

For e - 0, the steady state transfer functions from p and 0 l to 
w can be found. They are 


P 


7.2f + 0.9k 

3 * 6f 0 k 2 


a 


4f o + k i 

3 ‘ 6f 0 k 2 



The total control, therefore, is 



AAA 

c i + c^B + c CC + c W 

i a w 


c (7.21 +0.9^) - c cj (4W + 3.61.k 2 c u 

C . 1 + — 

1 3.6 f k 

o 2 


The value of e for a constant acceleration command (equal to i 
for a constant rate command) is therefore 



y7.2V0.9k;L) - c q< 4 V k l> + 3 * 6f Q k 2 C (jq 


3.6! kc. 
0 2 i 


For stability of the characteristic equation of (1), k 2 > 0. 

In general c_ < 0, and the three terms in the numerator are there- 
fore additive. Low k^ can be seen to decrease the error but it 
also decreases the damping of the characteristic equation of (1) and 
causes a higher control noise. 
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For other estimators, the error coefficients can be found in a 
similar way but since there are many nonzero states in the estimator, 
the expressions become complicated. 
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